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I.  INTRODUCTION 


Computer  programs  for  the  analysis  of  electromagnetic 
coupling  between  a  thin  straight  conducting  wire  and  an 
aperture  in  an  infinite  conducting  plane  are  briefly 
described  and  listed  in  this  report.  The  aperture  is  of 
arbitrary  shape  and  size.  The  wire  is  of  finite  length 
(with  or  without  loads)  or  of  infinite  length.  The 
excitation  is  either  a  plane  wave  incident  from  the  opposite 
side  of  the  wire  or  TEM  voltages  applied  on  the  wire.  The 
current  distributions  in  the  aperture  and  on  the  wire  are 
computed.  In  addition,  for  the  case  of  TEM  voltage 
excitaton,  the  power  transmitted  through  the  aperture  is 
computed.  For  the  case  of  plane  wave  excitation,  ~we 
evaluate^an  equivalent  circuit  of  the  aperture  for  the 
transmission  line  mode  on  an  infinitely  long  wire  or  an 
arbitrarily  loaded  wire.  It  is  assumed  that  the  reader  is 
familiar  with  [1],  where  the  general  theory  and  the  method 
of  computation  are  given. 


II.  COMPUTER  PROGRAM  DESCRIPTION 

There  are  two  main  programs:  MAIN1  and  MAIN2.  MAIN1 
is  for  the  case  of  an  infinitely  long  wire  or  an  arbitrarily 
loaded  wire.  MAIN2  is  for  the  case  of  an  unloaded  wire  of 
finite  length.  As  shown  in  Table  1,  each  main  program  calls 
subroutines  INDATA,  GEOM ,  CURDIR,  IMP,  MATRIX,  CSMINV,  MUL, 
MUL1,  and  GAUSS.  IMP  calls  subroutine  SICI  and  function 


subprogram  P.  MATRIX  calls  subroutines  TPARM,  SCAINT, 
VECINT,  and  TADJ.  SCAINT  calls  subroutine  INTGRL,  VECINT 
calls  subroutine  LININT,  and  CSMINV  calls  subroutine  DTRMNT. 
All  the  main  programs,  subroutines,  and  function  subprogram 
are  described  briefly  as  follows. 


Table  1.  A  list  of  main  programs,  subroutines,  and 
function . 


INDATA 


GEOM 


CURDIR 


MAIN1 


SICI 


MAIN2 


TPARM 


MATRIX  <  SCAINT— INTGRL 


VECINT— LININT 


TADJ 


CSMINV— DTRMNT 


MUL1 


GAUSS 


V 
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(A)  Main  Programs  MAIN  1  and  MAIN2 

For  program  MAIN1,  inputs  are  defined  by  DATA 
statements.  An  example  is 

DATA  NNODE/22/ , NEDGE/4 1 / , NFACE/20/ , NA/ 1 9/ 

DATA  D/0.25/,WL/0.4/,NB/20/,RB/0.001/ 

DATA  HPHI/ ( - 1 ,0. ) / , HTHETA/ ( 0 . ,0. )/, PHI/90 ./, THETA/0 . / 

DATA  ALAMDA/1./ 

DATA  V 1 /( 1 . ,0.)/,V2/(0. ,0.)/ 

DATA  G1/(0. 875,0. )/,  G2/(0. 667, 0. ) / 

Where 


NNODEs 

the  total  number  of 

patching  used  for  the 

nodes  of 

aperture . 

the 

triangular 

NEDGEs 

the  total 

patching . 

number  of 

edges  of 

the 

triangular 

NFACE= 

the  total 

tr iangul ar 

number  of 

patching . 

faces  (patches) 

of 

the 

NAs 

the  total 

number  of 

internal 

edges 

of 

the 

triangular  patching. 


m 


$1 


m 

% 

#1 


the  total  number  of  current  expansion  functions  on 
the  wire. 


the  distance  (in  wavelengths)  from  the  wire  to  the 


conducting  plane 


the  length  (in  wavelengths)  on  the  wire  where  the 
evanescent  current  is  assumed  to  exist. 


the  radius  (in  wavelengths)  of  the  wire. 


HPHI=  the  (^-component  of  incident  magnetic  field  (in 
units  of  ampere/meter) . 


HTHETAs  the  0-component  of  incident  magnetic  field  (in 
units  of  ampere/meter). 


PHI=  the  incidence  angle  <J>  (in  degrees) 


THETA=  the  incidence  angle  0Q  (in  degrees) 


ALAMDAs  the  wavelength  (in  meters). 


VI  ,V2= 


the  amplitudes  of  TEM  voltages  (in  units  of  2ZQ,  ZQ 
is  the  characteristic  impedance  of  the  trasmission 
line  formed  by  the  wire  and  the  conducting  plane) 
propagating  in  the  +z  and  -z  directions, 

respective! y . 
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G1  ,  G2=  the  reflection  coefficients  and  Tj*  (G1=  G2=  0. 

when  the  wire  is  matched  to  ZQ  or  is  infinitely 
long.  ) 


Outputs  of  program  MAIN1  are  the  coefficients  of 
current  expansion  functions  in  the  aperture  and  on  the  wire, 
the  power  transmitted  through  the  aperture,  and  network 
elements.  They  are  stored  in  a  data  file  with  I/O  unit 
number  21. 

The  minimum  allocations  are  given  by 

Complex  Y(NA,NA) ,Z(NB,NB) ,T(NA,NB) ,TT(NB,NA) ,YT(NA,NB) , 

GZ ( NB , NB ) , DGZ (NB,NB) , VI (NB ) ,VM(NA) ,ZZ(NB,NB) , 
YI(NB) , YTI (NA ) , CIA( NA ) ,CIB(NA) ,YY(NA,NA) 

Integer  NCONN(NEDGE, 3) , ITRAK(NEDGE) , IMIN(NEDGE) 


In  program  MAIN2,  inputs  are  the  same  as  those  in 
MAIN1.  However,  inputs  G1  and  G2  are  assigned  to  be  zero  in 
the  program,  and  NB  is  defined  as  the  sum  of  2  and  the  total 
number  of  current  expansion  functions  on  the  wire.  WL  is 
defined  as  the  difference  of  L  (the  length  of  the  wire)  and 
the  length  of  one  subsection,  i.e.,  WL=L/( 1+1 /NB) .  Ouputs 
are  the  same  as  those  in  MAIN1.  However,  the  TEM  equivalent 
circuit  is  not  evaluated.  The  minimum  allocations  are  given 
by 


ly.v/AV.v, 


Complex  Y(NA , NA ) , Z( NB , NB) ,T(NA,NB) , TT(NB , NA ) , YT(NA , NB-2 ) , 


GZ  ( NB-2 , NB-2 ) , DGZ(NB-2,NB-2) , VI (NB-2 ) ,  ZF( NB-2 , NB-2 ) , 
TTF(NB-2, NA) ,VM(NA) , ZZ( NB-2 , NB-2 ), YT (NA ) ,  YTI(NA) , 
CIA(NA) ,CIB(NA) , YY(NA , NA ) , TF( NA , NB-2 ) 

Integer  NC ONNCNEDGE, 3 ) , ITRAK (NEDGE ) , NBOUND(5 0 , 4 ) , IM IN ( NEDGE ) 

(B)  Subroutines  INDA TA ,  GEOM ,  and  CURDIR 


Subroutine  INDATA( DATNOD,  NCONN ,  NNODE ,  NEDGE)  reads 
two  sets  of  input  data  from  a  file  with  I/O  unit  number  20. 
The  subroutine  then  arranges  these  input  data  in  a  numerical 
order  for  the  triangular  patching  of  the  aperture.  The 
first  set  of  data  contains  node  numbers  along  with  their 
coordinates.  This  information  is  stored  in  an  output  array 
DATNOD.  The  second  set  contains  edge  numbers  with  the  node 
numbers  connected  by  them.  This  information  is  stored  in  an 
output  array  NCONN.  Note  that  in  the  input  data  file,  we 
enumerate  the  internal  edges  fist.  That  is,  if  there  are  N 
internal  edges,  the  edge  numbers  of  internal  edges  start  at 
1,  while  the  boundary  edges  start  at  N+1.  The  input 
variables  are  NNODE  and  NEDGE,  which  are  defined  in  (A). 
The  minimum  allocations  are  given  by 


Real  DATNOD(NNCDE,  3) 
Integer  NCONN (NEDGE, 3 ) 


Subroutine  GEOM (NCONN ,  NBOUND,  ITRAK,  IMIN,  NEDGE)  uses 
the  informations  stored  in  the  input  array  NCONN  to  form 


triangular  patches.  Inputs  are  NCONN  and  NEDGE.  Outputs 
are  array  NBOUND  storing  face  (patch)  numbers  and  their 
associated  edge  numbers.  ITRAK  and  IMIN  are  auxiliary 
arrays  needed  in  the  program.  Minimum  allocations  are  given 
by 

Integer  NC  ONN(NEDGE,  3),  NB0UND(50,  l»),  ITRAK(NEDGE )  , 
IMIN(NEDGE) 


Subroutine  CURDIR (NCONN ,  NBOUND,  NF  ACE,  NEDGE,  IMIN, 
NSE)  arranges  the  informations  in  input  arrays  NCONN  and 
NBOUND  and  then  transfers  them  into  an  output  data  file  with 
I/O  unit  number  21.  Inputs  are  NCONN,  NBOUND,  NFACE,  and 
NEDGE,  which  are  defined  previously.  Input  NSE  is  the  total 
number  of  boundary  edges.  IMIN  is  an  auxiliary  array. 
Minimum  allocations  are  given  by 

Integer  NCONN  (NEDGE,  3)»NB0UND(50,iO,  IMIN  (NEDGE) 


(C)  Suboutines  TPARM  and  TADJ 


Subroutine  TPARM(  N ,  DATNOD,  NCONN ,  NBOUND,  NNODE , 
NEDGE,  EN,  NN,  XN ,  ZN,  LN)  finds  the  edge  numbers,  node 
numbers,  the  x  and  the  z-coordinates  of  nodes,  and  lengths 
of  edges  for  the  triangle  whose  face  number  is  N.  These 
informations  are  stored  in  output  arrays  EN,  NN,  XN ,  ZN,  and 
LN,  respectively.  Inputs  are  N,  DATNOD,  NCONN,  NBOUND, 
NNODE,  and  NEDGE.  Minimum  allocations  are  given  by 

Integer  NCONN(NEDGE,  3) 

Real  DATNOD(NNODE, 3) 

Subroutine  TADJ(N ,  NA,  EN,  NN,  NCONN,  NEDGE,  M,  DIR) 
finds  the  current  reference  direction  crossing  an  edge  of  a 
triangle,  and  the  expansion  function  number  associated  with 
this  edge.  They  are  stored  in  output  variables  DIR  and  M, 
respectively.  DIR=1  if  the  current  direction  is  away  from 

the  triangle.  DIR=-1  if  the  current  direction  is  towards 

the  triangle.  DIR=0  if  the  edge  is  a  boundary  edge.  The 

edge  is  specified  by  N  which  is  the  node  number  of  its  free 

node  (node  not  on  the  edge).  Inputs  are  N,  NA,  EN,  NN, 
NCONN,  and  NEDGE.  Minimum  allocations  are  given  by 

Integer  NCONN(NEDGE, 3) 


(D)  Subroutine  MATRIX 


V//rf  r 
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Subroutine  MATRIX(NA ,  NB,  NNODE ,  NFACE,  DATNOD,  NCONN, 
NBOUND,  CIA,  CIB,  T,  TT,  Y)  computes  source  vectors  I*a  and 
tlb  and  matrices  [T],  [T],  and  [Y],  which  are  defined  in 
[1].  These  vectors  and  matrices  are  stored  in  output  arrays 
CIA,  CIB,  T,  TT,  and  Y,  respectively.  Inputs  are  NA,  NB, 
NNODE,  NEDGE,  NFACE,  DATNOD,  NCONN,  and  NBOUND. 


The  minimum  allocations  are  given  by 


Complex  CIA(NA) ,CIB(NA) ,T(NA, NB) ,TT(NB,NA) ,Y(NA,NA) 


Real  DA TNOD( NNODE , 3 ) 
Integer  NCONN (NEDGE, 3) 


(E)  Subroutines  IMP,  SICI,  and  Function  P 


Subroutine  IMP(NB,  RB,  Z)  computes  the  impedance  matrix 
[Z]  for  an  infinitely  long  wire.  Inputs  are  NB  and  RB,  and 
output  is  the  [ZJ  matrix  stored  in  an  array  Z.  Minimum 
allocations  are  given  by 


Complex  Z(NB,NB) 


Subroutine  SICI(SI, Cl, X)  computes  the  sine  and  cosine 
integrals 


SI  = 


sin  u 


Cl  r 


COS  u 


u 


du 


where  X  is  the  input,  and  SI  and  Cl  are  the  outputs 
subroutine  is  described  in  [2]. 


This 


Complex  function  P(AL,  Z,  ZL)  evaluates  the  scalar 
Green's  function 


P  = 


4  ttAI 


-jkJ(zm-z'  )2+  q- 


"  A1 


J(VZ,)2+  p2 


dz ' 


The  evaluation  is  described  in  [3].  Inputs  are  defined  as 


AL=  0.5Aln  (in  wavelengths). 

Z=  |zm”z'l  ^in  wavelengths),  the  distance  between  the 
field  point  and  the  midpoint  of  a  source  element. 

ZL=  p,  the  transverse  coordinate  (in  wavelengths)  of  the 
field  point. 


The  output  is  P. 


(F)  Subroutines  SCAINT,  VECINT,  LININT,  and  INTGRL 


1 1 

Subroutine  SCAINT(XS,  ZS,  X,  Z,  CPHI,  AREA)  computes 
the  integral  over  a  source  triangle  of  area  Aq  , 


CPHI 


where  R  is  the  distance  between  a  field  point  and  a  source 
point  in  the  triangle.  Inputs  are  the  x  and  the  z 
coordinates  of  nodes  of  the  triangle  and  of  the  field  point. 
These  coordinates  are  stored  in  arrays  XS  and  ZS,  and 
variables  X  and  Z,  respectively.  Aq  is  stored  in  the  input 
variable  AREA.  The  Output  is  CPHI. 


Subroutine  VECINT(XS,  ZS,  X,  Z,  CAXSI,  CAETA,  AREA) 
computes 


CAXSI 


CAETA 


where  5  and  n  are  area  coordinates  of  the  source  triangle 
defined  in  [1].  Inputs  are  XS,  ZS,  X,  Z,  and  AREA.  Outputs 
are  CAXSI  and  CAETA. 


MV 
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Subroutine  LININT(XS,  ZS,  X,  Z,  POTXSI,  POTETA,  AREA) 


computes 


POTXSI  = 


' 

POTETA  =  -£dS 
J  J  K 


Inputs  are  XS,  ZS,  X,  Z,  and  AREA,  and  outputs  are  POTXSI 
and  POTETA. 

Subroutine  INTGRL(XS,  ZS,  XF ,  ZF,  POT)  computes 


POT  = 


Inputs  are  XS,  ZS,  XF,  ZF,  and  AREA.  XF  and  ZF  are  the  x 
and  the  z  coordinates  of  the  field  point.  The  output  is 
POT. 


(G)  Subroutines  CSMINV,  DTRMNT,  MUL,  MUL 1 ,  and  GAUSS 


Subroutine  CSMINV(A,  NDIM,  N)  with  subroutine  DTRMNT 
inverts  a  NxN  matrix  stored  in  array  A.  The  result  is  also 
stored  in  array  A.  Inputs  are  A  and  NDIM=N. 


Subroutine  MUL(L,  M,  N,  A,  B,  C)  computes  the  product 
of  matrices  A(L,M)  and  B(M,N)  and  stores  the  result  in 
C(L, N) .  Inputs  are  integers  L,  M,  and  N  and  matrices  A  and 
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B.  The  output  is  C.  Minimum  allocations  are  given  by 
Complex  A(L,M) ,  B(M,N),  C(L,N) 


Subroutine  MUL1(L,  M,  A,  B,  C)  computes  the  product  of 
matrices  A(L,  M)  and  vector  B(M)  and  stores  the  result  in 
vector  C(L).  Inputs  are  integers  L  and  M,  matrix  A,  and 
Vector  B.  The  output  is  C.  Minimum  allocations  are  given 
by 

Complex  A(L,M)  ,  B(M)  ,  C(L) 


Subroutine  GAUSSCN,  A,  B,  EPS,  ISW)  solves  a  linear 
system  equation  AX=B  by  the  method  of  Gaussian  elimination. 
Inputs  are  the  matrix  A(N,  N),  N,  and  a  small  constant  EPS. 
The  output  ISW  =  1  if  the  absolute  value  of  the  pivot  of 
column  is  larger  than  EPS,  and  ISW  =  0  otherwise.  The 
solution  X  is  stored  in  B  as  an  output.  Minimum  allocations 
are  given  by 

Complex  A(N,N),B(N) 


uuuuuuuu 


III.  COMPUTER  PROGRAM  LISTING 
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ABOUT  Z-0,  OB  PRODUCING  SYMMETRIC  TEM  CUTHARD  TRAVELING 
CURRENTS)  FOR  THE  TEM  MODE  ON  THE  HIRE. 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

COMPLEX  Y ( 19 ,19) ,Z(62,62) ,T(19,62) ,TT(62,19) ,HPHI,HTHETA 
COMPLEX  YT ( 1 9 , 6 2 ) ,GZ (62, 62) ,DGZ  (62,62) ,?I (62) 

COMPLEX  VM  (19)  ,ZZ  (62,62)  ,YI  (19)  ,AI1  (2)  ,AI2{2)  ,Z  1  ,Z2,  YTI  (  19) 
COMPLEX  E5,V1,V2,JK,CIA  (19) ,CIB(19) ,G1,G2,PT,YY  (19, 19) 

COMPLEX  C1,C2,C3,C4,C5 

INTEGER  NCONN  (41,3)  ,ITHAK(41)  ,  NBOUND  (50,4)  ,IRIN(41) 

COMPLEX  VE1,VE2 
REAL  DATNOC (22,3) 

EQUIVALENCE (GZ,ZZ) 

COMHON/KKK/AK,  PI 
COMMON/JKK/JK 
COHMON/LHIRE/E,  WL 

COMMON/FIELD/HPHI,HTHETA,ERI, THETA 

COHMON/H  EU/K E, HU 

COHMCN/LCAD/G 1 ,G2 

COMMON/VOLT/ V 1, V  2 

COHMON/ZCO/ZO 

DATA  NNODE/22/, NEDGE/41/,NFACE/20/,NA/19/ 


DATA  D/0. 1/.81/1. 5/,NE/62/,RE/0.CG1/ 

EATA  HP H 1/  (-  1,0.)/,HTHETA/(0.,0.) /, PHI/90./, THETA/90./ 
DATA  ALAHDA/1.0/ 

DATA  VI/  (0,0)/,V2/<0.,0.)/ 

EATA  G1  / (0.8 75 ,0)/,G2/  (0.6667,0)  / 

FI=3. 14159265 
AK  =  2.*PI/ALAHEA 
JK= (0. , 1 . ) * AK 
VEL=3.E0£ 

AOMEGA=AK* VEI 
EPSLCN=1.E-09/(36.*PI) 

AHU=4.*PI*1.  E-07 
HE- ACM  EGA*EPSLON 
HU= AO  MEG  A*  A  HU 
US E=NEEGE-NA 
EH  1= PHI* FI/ 180. 

THETA=THETA*PI/ 1 80. 

Z0=60.*ALCG(2.*D/RB) 

OPEN  (UNIT=20,FILE=«IN.DATM 
OPEN  (ONIT=21,FIIE='COT. EAT») 

CALL  INDATA  (EATNOE,  NCO ti N  ,  N NCD E ,  NEDGE) 

CALL  GEOH(NCCNN,NBOaND,ITBAK, IMIN,NEDGE, NNODE) 

CALL  CUB  DIB  (NCONN, NBCUND,NPACE, HEDGE, ININ, NS E) 

HR ITE  (21,100)NA,NB,BB,HL,E 

FIND  IMPEDANCE  MATRIX  OF  HIR E,Z  (NB , NB)  ; 

CALL  IMP  (NB,EE,Z) 

FIND  COUPLING  MATRICES  T,TT, ADMITTANCE  Y, SOURCE  CIA,CIB 

CALL  MATRIX  (NA,NB, NNODE, NE DGE, NFACE, £ ATNOD, NCONN , 

1  NBOUN  D,CI A , CIB, T, TT , Y) 

DO  10  N=1,NA 
EO  10  N=1,NA 
YY(M,N)  =  Y(M,N) 

CALL  CSMINV  (Y,NA,NA) 

MATRIX  CALCULATIONS : 

CALL  MUL  (NA,NA,NB,Y,T,YT) 

CALL  MUL  |NB,NA,NB,TT,YT,GZ) 

EO  20  N=1,NB 
DO  20  N=1,NE 
EGZ (M,N) =  GZ  (  H,  N)  -Z  (  N,  N) 

K=1:TEM  INCIDENT  ;  K=2s  PLANE  HAVE  INCIDENT 

HBITE  (21,110) 

DO  30  K=1,2 
DO  32  M=1,NA 
VM (M) =-CIB  |H) 

IF  (K  -EQ.  1)  GO  TO  35 
HR ITE  (21 ,120) 

EO  34  M= 1, NA 
VH  (M)  =CI A  (M) 

CALL  MOL  1(NA,NA.Y,VN,YI) 

CALL  MUL1 (NB,NA,TT,YI,VI) 

EO  40  M= 1, NB 
DO  40  N» 1 , NB 


40  ZZ (H, N) =CGZ  (M,N)  16 

CALL  GAUSS  (HE, ZZ, VI, 11-11, ISM) 

IP  (ISM  .EQ.  1)  GO  TO  45 
TYPE  101 

101  FOB  HAT  ISW=Q  ,  STOP  * ) 

STOP 

45  BRITE  (21 ,130) 

NRITE(21,  140)  {VI  (H)  ,H=1,NE) 

AI  1  (K)=VI  (1) 

AI2(K)  =VI  (NB) 

CALL  HUL1  (NA,NE,YT,VI,YTI) 

CO  60  H=  1, NA 
VH  (H)  =  f  I  (N)-YTI(H) 

VH(H)=VH(H)/(120.*PI) 

60  CONTINUE 

BRITE(21, 150) 


BRITE  (21 ,140)  (VH  (H)  ,H=1,NA) 

BRITE (21,220) 

PINO  HIGH  NODE  ELE  CURRENT: 

Bl=(NB-2)/2 
DO  50  H=2,NU1 

ZN=BL*FLCAT ( 2*H-NB- 1) /FLOAT (2*NB-4) 

Hl=NB-fl*1 
E5=CEXP (JS*ZN) 

VI  (fl)=VI  (H)-VI  (1)  *E5 
VI  (HI)  =VI  (HI) —VI  (NB)  *P5 
CONTINUE 

BRITE  (21,140)  (VI  (H)  ,H=2,NB-1) 

FIND  TIHE-AVE.  POBER  TRANSfl.  THRO,  APERTURE, IT 

IF  (K.  EQ.  2)  GO  TO  30 
ET=  (0. ,0. ) 

DO  80  H-1,NA 
YTI (B) * (0. ,0.) 

DO  90  N=  1, NA 

YTI  (H)  =  YY  (H,N)  ♦Vfl(N)  *YTI  (H) 

PT=PT* VH  (H)  ♦CCNJG  (YTI  (H)  ) 

PT*0. 5*R  EAL (ET) 

BRITE (2 1, 240) PT 
CONTINUE 

FIND  EQUIVALENT  NETBORK : Z 1 ,Z 2, VE 1 , V E2  (NOBNALIZED  TO  ZO  ) 

II 1  OR  VO=  1) 

15=2.  MI 1(1)  *AI2  (1) 

21*-  (All  (1)  »lkI2  (1))/E5 

IF  (CABS  (All  (1)-AI2(1))/CABS  (All  (1)  )  .GT.  1.E-02)  GO  TO  52 
Z1=-2.*AI1  (1)*ZQ/(1.MI1  (1)) 

VE  1*-A1 1  (2)*(2.*Z0+Z1) 

BRITE  (2 1,230)  Z1,VEl 

GO  TO  17 

CONTINUE 

Z2=2.*(1.*AI2(1)  )/E5/(Al1(1)-AI2(1)) 

TO  AVOID  ANY  POSSIBLE  SHALL  ERROR  (SHALL  NEGATIVE  RESISTANT) 
DOE  TO  THE  VERY  SHALL  BCUNDCFF  ERROR  IN  REAL  (I  (  1) ) ,REA£ (I (NB) ) 


55  CONTINUE 

BZ  =  REAL(Z1)/CAES  (Z  1)  17 

IP  (  RZ  -GT.  0.)  GO  TO  25  ' 

IP(  ABS  (RZ)  -GT.  3.E-03  )  GC  TO  27 
Z1=(0.,1.)*AIHAG(Z1) 

25  BZ  =R  E  AL  (Z  2)  /CABS  (Z2) 

IF  (  RZ  -GT-  0.  ) GO  TO  37 

IF  (ABS  (RZ)  -GT.  3.E-C3)  GC  10  27 

Z2=(0-,1-)*AIMAG(Z2) 

37  CONTINUE 

E5=1.*Z1+Z2 

VE1=-E5*AI 1(2) *Z2*AI2 (2) 

VE2=E5*AI2  (2) -Z2*A1 1  (2) 

WHITE  (21, 170)  Zl,Z2, VE1.VE2 
GO  TO  17 

27  WRITE (21,2  10) 

100  FOBMAT(//,2X,*NA=*,I4, 1X,*NB=* ,14, 1X,*RB=* ,F6.4,  IX, 

1  *WL=*,F6.4, 1X,*D=*,F6.4) 

110  FORMAT  (/,2X,«TEM  WAVE  INCIDENT*) 

120  FORMAT (/,2X, ‘PLANE  WAVE  INCIDENT*) 

130  FORHAT(2X,*EIEC.  CURRENT  ON  WIRE, 1ST  6  LAST  ARE  TEH  AT  Z~0-,0»*) 
140  FOBMAT (2X,2E) 

150  FORHAT(2X,*MAG.  CURRENT  DENSITY  ON  APEBTUBE/ (120. *PI) =*) 

170  FOBMAT (/, 2 X, *  EQUIVALENT  NETWORK (NORMALIZED  TO  ZO) :  ' ,/, 2X, • Z 1= • , 
1  2E,2X,*Z2=*,2E,/,2X,» VE1=* ,2E,2X,*VE2=* ,2E) 

210  FORMAT  (2X,**«*  STOP,  NEED  LARGER  WL  S  NB  ***') 

220  FORHAT(2X,* HIGH  MODE  ELEC.  CURRENT  ON  WIRE*) 

230  FOBM  AT (2X,  * Z2  IS  OPEN*,/,*  Z  1=  •  , 2E, 2X, *  VF  1=*  2E) 

240  FORMAT  (2X, 'POWER  TRANSH.  THRO.  APERTDHE,PT=* ,E) 

17  CLOSE  (UNIT=20,FILE=*IN.CAT») 

CLOSE  (UNIT=21,FILE5=*CUT.DAI*) 

STOP 

END 


ftftnnftnftflftftnnnnoftnoftftftftft 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
BA  IN  2: MAIN  PROGBAN  FOB  TEE  PBOBLEH  CF  AN  UNLOADED  VIBE  OF 
FINITE  LENGTH  BEHIND  AN  AFEFTUBE  OF  ABBITBARI  SHAPE  AND  SIZE. 
THE  APERTURE  IS  IN  AN  INFINITE  CONDUCTING  PLANE 
OF  ZERO  THICKNESS. 

EXCITATION  IS  EITHER  A  PLANE  VAVE  FRCH  THE  OEPCSITE  SIDE 
OF  THE  HIRE  OR  TEH  VOLTAGES  ON  THE  VIBE. 

IN  PUTS :N NODE ,N EDGE, NFACE, NA=THE  TOTAL  NUHEEBS  CF 
NODES, EDGES, FACES  (PATCHES) , INTERNAL  EDGES 
OF  TRIANGULAR  PATCHING  FOR  IBE  PAERTUBE. 

NA  IS  ALSO  THE  TOTAL  NUHBER  OF  EXPANSION 
FUNCTIONS  IN  THE  APERTURE. 

D: DISTANCE  BETWEEN  THE  HIRE  AND  THE  CONDUCTING  PLANE 
VL: LENGTH  OP  MIRE  -  LENGTH  CF  ONE  SUBSECTION 
NE:TOTAL  NUMBER  OF  EXPANSIONS  ON  THE  HIRE  ♦  2. 

HPHI,HT BETA, PHI, THETA: THE  INCIDENT  H-HILEDS  AND  ANGLES. 
ALANDA:NAVELENGTH  (IN  HETERS) 

VI, V2: AMPLITUDES  OF  TEH  VOLTAGES  (IN  UNITS  OF  ZO) 

RB: RADI  US  OF  THE  HIRE 

OUTPUTS:  COEFFICIENTS  OF  CURRENT  EXPANSION  FUNCTIONS 
IN  THE  APERTURE  FOB  BOTH  EXCITATIONS. 

CURBENTS  ON  THE  HIRE 

THE  TIME-AVE.  POWER  TRANSMITTED  THROUGH  TH  APERTURE, 

FOB  TEN  EXCITATION. 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

COMPLEX  Y ( 19, 19) ,Z( 11, 11) ,1(19,11) ,11(11,19) ,HPHI,HTHETA 
COMPLEX  YT  (19,9)  ,GZ{9,9)  ,DGZ  (9,9)  ,VI(9) 

COMPLEX  ZF(9,9),TF(19,9)  ,TTF(9,19) 

COMPLEX  VH  ( 1 9)  ,ZZ  (9,9)  ,YI  (19)  ,  All  (2)  ,  AI2  (2)  , Z  1,Z2, YTI  (  19) 
COMPLEX  E5,V1,V2,JK,CIA(19),CIB(19) ,G1 ,G2 ,PT,IY (19 , 19) 

INTEGER  NCONN  (41,3)  ,  ITB  AK  (41)  ,  NBOUN  D  (50 , 4)  ,IHIN(41) 

BEAL  DATNOD  (22,3) 

EQUIVALENCE (GZ,ZZ) 

COMMON/KKK/AK ,PI 
COMHON/JKK/JK 
COMHON/LHIRE/D, VL 

COHMON/FIEL D/HPH I , HTHETA ,P HI,  THETA 

COMHON/H EU/V E , HU 

COHHON/VOLT/Vl,V2 

COHMON/ZOO/ZC 

COMMON/LCAD/G 1 ,G  2 

DATA  NNODE/22/, N EDGE/4 1/,NF AC E/20/, NA/ 19/ 

DATA  D/0.25/, HL/ 0.45/, NB/1 1/,RB/0.00 1/ 

DATA  HPHI/ (- 1, 0. )/, HTHETA/ (0.,0.)/, PHI/90./, THETA/90./ 

DATA  ALAMDA/1.0/ 

DATA  Vl/(1,0.)/,V2/(0.,0.)/ 

G1  =  (0,0) 

G2=  (0,0) 

KBF=NB-2 

FI=3. 14159265 

AK=2.*PI/ALAMDA 

JK=  (0. ,  1.) *AK 

VEL-3. E08 

AOMEGA=AK*VEL 

EPSLON=1. E-09/ (36. *  PI) 

AHU*4.*PI*1.E-07 

SE»AOMEGA*EFSLCN 

VU*AOMEGA*AMU 
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BSE=NEDGF-NA 

FBI=PHI*FI/180. 

THETA=THETA*EI/180. 

OPEN  (UNIT— 20,  FILE^’IN.DAT') 

OPEN  (0NIT=21,FILE='CUT.EAT') 

CALL  IN  DATA (EAT NOD, NCONN, NNCDE, NEOGE) 

CALL  GEOB  (NCO NN,  NBOU N D,ITRAK, IMIN, N EDGE,  N NODE) 

CALL  COR DIB  (NCONN, NBCUND, NFACE, NEDGE ,IHI N , NSE) 

1BITE  (21, 100) NA.NB,BB,ML,E 
C 

C  FIND  IflPIDANCE  MATRIX  OF  NIRE,Z  (NB, ME)  ; 

C 

CALL  IMP  (NB,BE,Z) 

C 

C  PIND  COUPLING  MATRICES  T ,TT , A DMITT A  NCI  I, SOURCE  CIA,CIB 

C 

CALL  MATRIX  (N A, NB,NNOCE, NIDGE, NFACE, E1TNOE, NCONN, 

1  ABOUND, CIA, CIB,T,TT,Y) 

EO  10  fl=  1 ,  NA 
DO  10  N=  1,  NA 
10  IT(H,N)=Y(M,N) 

CALL  CSMINV (Y,NA,NA) 

C 

C  BEDUCE  Z,T,TT  TO  ZF,TE,TTF  FOR  UNLOADED  WIRE: 

C 

DO  500  M=  1  ,NEF 
DO  510  N= 1 ,NEF 
510  ZF(M,N)=Z(N+1,N*1) 

DO  500  N= 1 , NA 
TTP(M.N)  =TT(B*1,N) 

500  TF(N.N)=-TTF  (M,N) 

BB=NBF 

C 

C  MATRIX  CALCULATIONS: 

C 

CALL  MUL  (NA,NA,NE,T,TF,YT) 

CALL  MUL  (NE, NA,NB,TTF,YT,GZ) 

EO  20  H=  1*  NB 
DO  20  N=  1 ,  NB 

20  DGZ(M.N)  =GZ(fl,N)-ZF(N,N) 

C 

C  ' K= 1 :TEM  INCIDENT  ;  K=2:  PLANE  NAVE  INCIDENT 

C 

Z0=60.*ALOG ( 2. *D/R  B) 

VBITE  (21 , 1 10) 

DO  30  K=  1,2 
DO  32  M= 1 , NA 
32  ?IKH)=-CIB(M) 

IF  (K  .EQ.  1)  GO  TO  35 
MBIT!  (21,120) 

DO  34  11=1, NA 

34  VH (fl) =CIA  (M) 

35  CALL  M0L1  (NA,NA,T,VM,YI) 

CALL  MUL  1(NB,NA,TTF,YI,VI) 

DO  40  fl=1,NB 

DO  40  N= 1, NB 
40  ZZ  (H,M)=DGZ  (R,N) 

CALL  GAUSS(NE,ZZ, VI, IE-11, ISB) 

IF  (ISN  .EO.  1)  SO  TO  45 


101 


IS  W=0 ,  STOP  * ) 


20 


45 


60 


C 

c 

c 


90 

80 


30 

100 

110 

120 

130 

140 

150 

240 

17 


FORMATS 
STOP 

WRITE  (21.130) 

WRITE  (21. 140)  (VI(8)  , (1=1. RE) 
CALL  BUL1(NA,NE,YT,VI,YT1) 
DO  60  8=1,  NA 
VB  (8)  =YI  (M)-YTI  (8) 

VB  (M)=VM(B)/(120.*PI) 

CONTINUE 

WRITE (21, 150) 

WRITE  (21 ,140)  (Vfl  (M)  ,  8=1 ,  NA) 


FIND  TI8E-AVE.  POWER  TRANSH.  THRO.  APEBTURE.PT 

IF (K. EO. 2) 60  TO  30 
PT= (0- , 0- ) 

DO  80  8=  1 ,  NA 
YTI  (8)  =  (0.  ,0.) 

CO  90  M=  1 , NA 

YTI  (8)  =Y  Y  (8,  N)  *V8(N)  *YTI  (8) 

PT=PT*VB (8) *CCNJG (YTI (8)  ) 

PT=0. 5*REAL (FT) 

WRITE (21 ,240) PT 
CONTINUE 

F0E8AT  (//,2X,*NA—  *,I4,1X,'NB=',I4, 1X,'HB=' ,F6.4,1X, 

1  •WL=* , F6.4,1X,'D=',F6.4) 

FORMAT  </.2X, *TEM  WAVE  INCIDENT') 

FORMAT  (/,2X, 'PLANE  WAVE  INCIDENT') 

FOFBAT(2X,'ELEC.  CURRENT  CN  WIBE') 

FORMAT  (2X,2E) 

FORMAT  (2X, 'BAG.  CURRENT  DENSITY  ON  APERTURE/ (120.*PI) =• ) 
FORMAT (2X, 'POWER  TRANSB.  THRO.  APERTURE ,PT=* .E) 

CLOSE  (ONIT=20.FILE='IK.DAI') 

CLOSE  (DNIT=2 1 , FILE=' OUT. E  AT* ) 

STOP 

END 


21 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
C  INDATAtREADS  two  sets  of  input  data  pcb  the  c 

C  TRIANGULAR  PATCHING  OF  THE  PAERTOBE: ( 1)  NOCE  NOS.  C 

C  WITH  COORDINATES,  STORED  IN  DAI  NOD.  (2)  EDGE  NCS.  C 

C  AND  NODE  NOS.  CONNECTED  EY  THEN,  SOTRED  IN  NCONN  C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  I NCAT A  ( D  ATNCD  ,  NCCNN  ,  N BCD E  ,  N  E  DGE) 

DIMENSION  DATNOD(NNODE,3) 

INTEGER  NCONN  (NEDGE, 3) 

DO  10  I=1,NNCDE 
READ  (20,5) NOCE, X, 2 
5  FORMAT  (12, 2F9.6) 

AN  =  FLOAT  (NCCE) 

CATNOD  (NODE,  1)  =  AN 
DATNOD  (NCDE,2) =  X 
CATNOD  (NCDE,  3)  =2 
10  CONTINUE 

DO  20  1=  1, NEEGE 
BEAD  (20,  15)  NE  ,  NF,  NT 
15  F0RMAT(I3,2I2) 

NCONN  (NE,  1)  =NE 
NCONN  (NE,  2)  =NF 
NCONN  (NE,  3)  =  NT 
20  CONTINUE 

8RITE  (2 1 ,  18) 

18  FORM  AT  (*  1*  ) 

WRITE  (21,  19) 

19  FORMAT  (6X,*  TEE  FOLLOWING  IS  THE  INFORMATION  CONCERNING  NODES 
1  AND  THEIR  COORDINATES* ,/) 


DO  30  1= 1, NNCDE 
IDUMMY=IFIX (DATNOD (1,1)) 

WRITE  (21,  21)  IDUMMY, DATNOD  (1,2)  , DATNOD  (1,3) 

21  FORMAT  (3X,» NOCE  NUMBER-* ,13 ,3X , * X-COOFCINATE-* ,F7.4,3X, 

1*  Z—COORCINATE-*  ,F7.4) 

30  CONTINUE 

HR  IT  E  (2 1 ,  28) 

28  FORMAT  (*1*) 

WRITE (21 ,  29) 

29  FORM  AT (1  OX, *  THIS  IS  THE  INFORMATION  CONCERNING  EDGES  ,  NODES* , 
DO  40  1=  1, NEDGE 

WRITE (21,  31) NCONN (1,1) , NCONN (I, 2) , NCCNN (1,3) 

31  FORMAT(3X,*EDGE*,I3,1X,*IS  CONNECTED  FROM  NODE* , IX, 13,1 X, 

1*  TO  NODE*, IX, 13) 

40  CONTINUE 

RETURN 
END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  GEOM: USES  INFORMATIONS  IN  NCCNN  TC  FORM  TRIANGULAR  C 

C  PATCHING. DIMENSION  OF  NECUND  MUST  BE  INCREASED  FCR  C 

C  TOTAL  NUMBERS  OF  FACES  >  50.  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  GFOM ( NCONN, NBCUND, ITRAK ,IMIN , HEDGE) 

INTEGER  NCONN  (NEDGE, 3) , NEOUND (50, 4)  ,ITRAK (NEDGE) 

INTEGER  IMIN  (NEDGE) 

COM MON/IF/I FACE 
IFACE=0 
NF  1=0 
NF  2=0 

DO  100  XJ*1, NEDGE 
3COONT=0 


VV:.*. 


i1 .  •  .  - 


N1=NC0NN  (I J , 2) 

82= NCONN  (I J,  3) 

CO  10  1=  1, NEDGE 
DO  10  J=2, 3 
IP (I.EQ. IJ) GO  TO  10 
8A  =  NCONN |I,J) 

IF (NA. EQ. N 1. OB.NA. EQ. N2) GO  TO  6 

GO  TO  10 

ICOUNT=ICOOMT ♦ 1 

ITBAK  (ICGUNT)  =1 

CONTINUE 

RABK 1=0 

HABK  2-0 

CONTINUE 

K 1=  1 

II =ITB AK f K 1 ) 

DO  15  1=2, ICOONT 
IF  (ITBAK <I). IT. 1 1 ) GO  TO  12 
GO  TO  15 
II =ITB AK (I) 

K 1  =1 

CONTINUE 

IF  (HARK1.EQ.  ICOUNT)  GO  TO  100 
IF  (I1.GT.IJ) GO  TO  20 
GO  TO  31 
CO NTINUE 
N3=NCONN (I 1, 2) 

84=NCONN  (11,3) 

IF(N3.EQ.Nl.CB.N3.EQ.N2)GC  TO  21 

1F(N4.EC.N1.CB.N4.EC*N2)GC  TC  22 

NB  =  N4 

GO  TO  23 

*B=N3 

CONTINUE 

ICO=0 

DO  25  1=  1 , NE EGE 
CO  25  J=  2,  3 
IF  (I.EQ.II)GC  TO  25 
8C=NCONN  (I,J) 

IF  (NC.EQ.NB)  GO  TO  24 

GO  TO  25 

ICG=ICO»1 

ININ (ICO) =1 

CONTINUE 

DO  30  1=1,  ICC 

IA=IHIN  (I) 

IF (Nl. EQ. NCONN ( I A  ,  2) .OB.N1.EQ. NCONN (I A , 3) ) GO  TC  29 

IF  (N2.  EQ.  NCONN  (I A ,2) .OB. N2. EC« NCONN (IA,3) ) GO  TO  29 

GO  TO  30 

12=1  A 

GO  TO  32 

CONTINUE 

COMTINUE 

ITBAK  (K1)  =  NE EG E*1 
8ABK 1=NABK 1*1 
GO  TO  75 

IF  (I2.LT.  IJ)  GO  TO  74 
IP  (IFACE.EQ.  0)GO  TO  33 
NF1=NB0UND(IFACE,2) 

NF2=NBOUND(IFACE,3) 


nnoooononoonn 
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33  IF (IJ.EQ.NF1.AND. 12. EQ. NF2) GC  TC  74 

IF ACE= IF  AC  E*  1 
NBOUND(IFACE,1)=IFACE 
NBOUNDfIFACE, 2)  =  I J 
ABOUND  (IFACE,3)  =  I1 
NBCU  ND ( IF  ACE,  4)  =12 
MARK2  =  MARK2«-1 
MARK1=MARK 1*1 


IF (MARK2. EC- 2) GC  TO  100 
ITRAK(K1)=N£DGE*1 
GO  TO  75 
74  CONTINUE 

ITRAK (K1)=NEDGE*1 
MARK 1=MARK 1*1 
GO  TO  75 

100  CONTINUE 

iR ITE  (21#  98) 

98  FORMAT ( *  1  * ) 

WRIT E  (2 1 #  99) 

99  FORMAT(10X,  'THIS  IS  THE  INFOR MATIC N  CONCERNING  PACES,  EDGES*,/) 
NR  IT  E  (2  1  #  101)  i (NBOUND(I,J) , J=1 ,4) ,I=1,IFACE) 

101  FORMAT  (3X, 'FACE*  ,13,  IX,  VIS  EETNEEN  THE  EDGES *,  1  X ,13 , 

11X,I3#1 X ,13 ) 

EO  120  1=  1 , IF  ACE 
EO  120  J=2,4 
ISEDGE=NEOUNE (I, J) 

KCOUNT=0 

EO  125  K=  1, IFACE 
EO  125  M=2,4 

IF (I.EO.K. ANC.J.EQ.M) GO  TC  125 
IF  (ISEDGE. EO. NEOUND  <K,K) ) NCOUNT=NCOUNT* 1 
125  CONTINUE 

C  IF  (NCOUNT-EQ-0)  HRITE  (21,  150)  ISEDGE 

120  CONTINUE 

150  FORMAT (/5X,I3,2X,' IS  A  BCUNEABY  ECGF*) 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


AIL  NODES  OP 


CP 


TPABM:OBTAINS  PARAMETERS  CF  A  TRIANGLE 
IN  PUT: N : TRIANGLE  NUMBER  CONSIDERED 

DAT  NOD  (  NNODE, 3) :NOEE  #,X.Z.COOE  FOR 
SYSTEM 

NCONN  (NNODE, 3) :ECGE  NO., NODE  NO.  (FRON-TO) 

N BOUND  (50,4) :FACE  f , EDGE  *S  FOB  ALL  TRIANGLES 
SYSTEM 

OUTPUT:EN(3) :ECGE  #  CF  THE  TRIANGLE  (SANE  ORDER  AS 
INPUT  EY  USER) 

NN(3)  :  NODE#  (ORDERED  AS  NODE  BETNN. EN  ( 1)  ,  (2)  ;  (2)  , 

(3)  ;  (3)  ,  (1) 

XN  (3)  ,ZN  (3)  :  X,  Z-  COCR.  OF  VEB1ICES 
L N  (3)  :  L  ENGT H  OF  EDGES  OPPOSITE  NODE  NN  ( 1)  ,  (2)  ,  ( 3) 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  TEARM  (N, EATNCC, NCCNN, NECUNC, NNOCE, N! CGI, 

1  EN ,  NN , XN , Z N, L N) 

INTEGER  EN  (3)  ,  NN  (3) 

INTEGER  NCONN(NEDGE,3) , NBC UND (50,4) 

REAL  LN  (3)  ,XN(3)  ,ZN(3)  ,  DAT  NOD  (NNOCE,  3) 

CO  10  M=  1, 3 

10  EN(N)=NBCUND(N,M*1) 

CO  20  M=  1, 3 


mt*  I.ViAifc.  Ii 


Wl 

m 


I1=EN(fl) 

I2=M*1  2U 

IF  (12. GT. 3)  12=1 
I2  =  EN (12) 

DM  (ii) =  NCON  N  (11,2) 

IF (NCONN (II, 2)  .EQ.NC0NN(I2,2))  GO  TC  20 
IF  (NCONN  (II, 2). EQ. NCONN  (12,3)  )  GO  TO  20 
NN(fl)=NCONN(I1,3) 

20  CONTINUE 

CO  30  H=1,3 

IN  (M)  =OATNCD  (NN  (M)  ,2) 

30  ZN  (M)  = DAT NOD  (NN  (fl)  ,  3) 

IN (1)  =  SQBT  (  (XN  (2)  -XN{3) )  **2*  (ZN(2)  -ZN(3)  )  **2) 

LN  (2)=SQRT{  (XN  (3) -XN  ( 1)  )  **2*  (ZN  ( 3) -ZN  ( 1))  **2) 

LN  (3)  =  SQBT  (  (XN  ( 1 )  —X N  (2)  )  **2  ♦  (ZN  ( 1) -ZN  (2) )  **2) 

EETURN 

END 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  TADJ: ADJUST  EDGE  NO.  S.T.  H  OPPOSITE  NOCE  N;  ADJUST  INDEX  H, 

C  S.T.  M  STARTS  AT  1  AND  EXCLUDES  BCUNDABI  EDGE;  OBTAIN  DIB,  THE 

C  DIRECTION  OF  THE  CURRENT  THRO  THE  EDGE,  CIB=1  IF  COBBENT  AHA Y 

C  TI1E  TRIANGLE, -1  IF  TOWARDS  THE  TRIANGLE, 0  IF  IT  IS  BCUNDAEY. 

C  INPUT: N-EOINTER  OF  NODE  NC  CONSIDERED  IN  A  TRIANGLE8 1 , 2, 3 

C  KSE:#  OF  SURFACE  (BOUNCABY)  EDGES 

C  EN  (3)  :  EDGE  *  OF  THE  TRIANGLE  (=  1 , 2,.  . ,  HEDGE) 

C  NN  (3) : NODE  *  OF  THE  TBIANGLE  , SANE  CBDEB  AS  IB  TPABH 

C  ITBAK(NEDGE): AUXILARY  VECTOR 

C  NCONN (NECGE, 3) : NODE 

C  NECGE: I  CF  EDGES 

C  OUTPUT: 

C  H: INDEX 

C  DIR: DIRECTION  OF  CURRENT  (=1,-1,0) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  TADJ(N,NA,EN,KK, NCO KN,NEDGE , H , DIR) 

INTEGER  NCONN  (NECGE, 3)  ,EN<3)  ,NN(3) 

EIB=0. 


C  ADJUST  EDGE  «,S.T.  H  DENOTES  THE  DEGE  «  OPPOSITE  NODE  N 

C  IF  H  IS  A  BOUNCABY  ECGE,J0HF  OUT 

C 

IF  (N  .EQ.  1)  M=EN  (3) 

IF  (N  .BQ.  2)  H  =  EN(1) 

IF  (N  -EQ.  3)  H=EN  (2) 

IF  (H.GT. NA)  RETURN 
C 

C  FIND  DIR: 

C 

N2=N+1 
N3  =  N+2 

IF (N 2  .GT.  3)  N2=N2-3 
IF  (N3  .GT.  3)  N3=N3-3 
cia=  i.o 

IF  (NCONN  (H,3)  .  EQ.  NN  (N2)  -AND.  NCONN  (H,2)  .  EQ.  NN(N3)  )  DIH=-  1. 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  VECINT: THIS, SIT H  LININT, EVA1UTES  VECTOR  POTENTIAL  INTEGBAl 

C  OVER  A  TBIANGLE  BEG10N 

C  IN PUT:COOR.  OF  3  VERTICES  OF  THE  TBIANGLE  (XS  (3)  ,ZS  (3)  ) 

C  CB VEBATION  POINT  (X,Z) ,  AREA  (CF  THE  TRIANGLE) 


N2=N2-3 

N3=N3-3 


|$ftSwK5< 


OUTPUT: INTEGRAL 


CAXSI ,CAETA 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SOBROUTINE  VECINT(XS,ZS,X,Z, CAXSI, CAE1A, ABEA) 

IMPLICIT  CONFLEX  (C) 

BEAL  CABS,  COS,  XS  (3)  , ZS (3) 

COMMON/KKK/AK,PI 
COMNON/V  EC/XS I  (7)  ,ETA(7) 

CF=CMPLX(0.,0.) 

CG-CNPLX  <0. , 0. ) 

DO  120  1=1,7 

B1  =  (  (X-XS  (  1)  )-(XS(2)-XS(1))  *XSI(I)-(XS(3)-XS{1))*ETA(I))4*2 
B2=  (  (Z-ZS  (1))-(ZS(2)-ZS(1))  ♦XSI  (I)  -  (ZS  (3)  -  ZS  { 1)  )  *ETA  (I)  )  **2 
B=SQRT  (R  1+R2) 

CB=Cf1PLX  {0.  0,-1. 0*AK*R) 

IF (CABS (CB )  .LE.  1.0E-06)GO  TO  102 
CA=  (CEXP(CR)  —CMP  LX  (1.0, 0.0))  /CMPLX  (B,0.) 

CF  1=CMPLX  (XSI  (I)  ,0.)  *CA 
CG1=CMPLX (ETI (I) ,0.  )  *CA 
GO  TO  103 

102  CF1=CHPLX(0. ,-AK*XSI (I)) 

CG  1=CMPLX(0.  ,—  A K * ET A  (I)  ) 

103  IF (I. EQ. 1 ) GO  TO  105 

IF  ( I.  EQ.  2  .OB.  I.EQ.3  .OB.  I.EQ.4)GC  TO  110 
CF=CF+CF1*CHPLX (0. 1259392,0.) 

CG=CG*CG1*CHPLX(0. 1259392,0.) 

GO  TO  120 

105  CF=CF*CF1*CMPLX (0.225,0.) 

CG=CG*CG1*CHILX  (0.225,0.) 

GO  TO  120 

110  CF=CF*CF 1*CH  PL X ( 0. 1323942, 0. ) 

CG=CG*CG1*CMPLX(0. 1323942,0.) 

120  CONTINUE 

CALL  LININT (XS,ZS,X,Z , POT XSI , POTET A, ABEA) 

CAXSI=CF<*CMPLX  (AREA,0. ) ♦CMFLX (POTXSI,0.) 
CABTA=CG*CMPLX(AREA,0.) ♦CMPLX (PCTETA,0.) 

150  CONTINUE 

BETUBN 
END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  SCAINTrTHIS  SUBROUTINE, HITH  SUBBCUTINE  INTGRI, EVALUATE 

C  THE  SCALB  POTENTIAL  INTEGRAL  OVEB  A  TIANGLE  BEGICN. 

C  IN PUT:COCB.  OF  VERTICES  CE  TBIANGLE  XS(3),2S(3) 

C  ABEA  OF  THE  TBIANGLE,  AREA; CBSEBVATION  POINT  (X,Z) 

C  AN, PI  IN  COMMON  KKK 

C  OUTPUT :SC ALAR  POTENTIAL  INTEGRAL  OVEB  ABEA, STORED  IN  CPHI 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  SCA INT (XS, ZS , X,Z ,CPHI, ABE A) 

IMPLICIT  COMPLEX  (C) 

BEAL  CABS,  COS,  XS(3)  ,  ZS  (3) 

COHMGN/K RK/AK ,PI 
COHHON/VFC/XSI (7) , ETA (7) 

XSI(1)  =  1. 0/3.0 
XSI|2) =0.05971587 
XSI(3)=0. 47014206 
XSI  (4)  =XSI  (3) 

XSI  |5)=0. 79742699 
XSI(6)=0.  10128651 
XSI  (7) =XSI  (6) 

ITA  (  1)  =XSI  (1) 

ETA  (2)  =XSI  (3) 


OTPO 
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ETA  (3)  =XSI  (2) 

ETA  (4)  =XSI  (4) 

ETA  (5)  =XS I  (6) 

ETA(6)=XSI  (5) 

ETA (7) =XSI  (7) 

CF=CHPLX  (0.0,0.  0) 

EO  120  1=1,7 

51=  (  (X-XS  (  1)  )  -  (XS  (2)  -XS  (  1)  )  *XSI  (I)  -  (XS  (3)  -XS  ( 1)  )  *ETA  (I)  )  **2 
B2=  ((Z-ZS(1)  )-  (ZS  (2)  -ZS  ( 1)  )  *XSI  (I)-(ZS(3)  -ZS  ( 1) )  *ETA  (I) )  **2 
B=SQBT(R  1  +  R2) 

CR=CHPLX  <0.0,-1.0*AK*B) 

IF(CABS(CR)  .IE.  1.0E-06)  GO  TO  102 

CF1=  (CEXP  (CR) -CHPLX  (1. 0, 0.  0)  ) /CNPLX (R,0. 0) 

GO  TO  103 

102  CF 1=CNPLX  (0.0,-AK) 

103  IF  (I  -EQ.  1)  GO  TO  105 

IF  (I.  EO.  2  .05.  I.  EQ.  3  .05.  I-EC-4)GO  TO  110 
CF=CF*CF1*CMPLX (0. 1259392,0.0) 

GO  TO  120 

105  CF=CF*CF1*CMPLX (0.225,0.0) 

GO  TO  120 

110  CF=CF+CF1*CMPLX(0. 1323942,0.0) 

120  CONTINOE 

CALL  INTGRL (XS,ZS,X,Z, POT) 

CPHI=CF*CHPLX (AREA, 0.0) +CNFIX (POT, 0.0) 

150  CONTINUE 

RETURN 
END 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  LININT: THIS, WITH  SUBBOOTINB  INTGBL ,EVAVUIATES  XSI/R  AND 

C  ETA/R  INTEGRALS  OVER  A  TRIANGLE  REGION. THE  QUANTITIES  DEFINED 

C  HEBE  ARE  THE  SANE  AS  THOSE  IN  THE  REFERENCE. 

C  AND  POT  1  FROM  INTGRI  IS  USEE. 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  LININT  (XS, ZS,X, Z , POTXSI, PGTET A, AR FA) 

BEAL  XS  ( 3)  ,ZS  (3) 

COHHON/POTEN/POT 1 

A=  <XS(2)  — XS(1))**2+(ZS(2)-ZS(1))**2 
8=  (XS(3)-XS(1))**2  +  (ZS  (3)-ZS  (1) )  **2 

C=-2.*((X-XS(1))*  (XS(2)-XS(I))  *(Z-ZS(1))  *  (ZS  (2) -ZS  (  1)  ) ) 

C=—2.*  ( ( X—  XS  fl) )  ♦  (XS  (3)-XS  (1)  )  ♦  (Z-ZS  (1) )  ♦  (ZS  (3)-ZS  ( 1)  )  ) 
E=2.*|(XS(2)-XS(1))*  (XS(3)-XS(1))  +  (ZS  (2) -ZS  (  1)  )  ♦  (ZS  (3)  -ZS<1)  )) 
F=(X-XS(1))**2+(Z-ZS(1))**2 

A1=  (2.*B-C+D-E)*SQRT(B+D  +  F)  ♦  (2.  *  A+C-D-E)  ♦  SQRT(A*C+F) 

A2=4.  *  (A+B-E) 

A3=A 1/A2 

A4=4.*  (A+C)  *  (B  +  D+F)  *4.*F*  (B-C-EJ-  (C  +  D+E)  **2 
A5=8.*SQRT( (A+B-E) **3) 

A6=A4/A5 

IF(ABS(A6).LE.1.  E-04)  GC  TO  5 
AL  1=2.*3QRT(A  +  B-E)  *SQRT  (E+D+F) 

AL2=2.*SQRT(A  +  E-E)*SCRT  (A+C  +  F) 

AL3=  (2.*E-C  +  E-E) 

AL4=  (2.*A+C-D-E) 

AJ  1=A3+A€*ALCG(ABS((AL1+AL3)/(AL2-AI4))) 
AJ3=A3+A6*ALCG(ABS((AL2  +  AL4)/(AL1-A13)) ) 

GO  TO  6 

5  AJ  1= A  3 


AJ3=A3 

El *S0RT (A+C+F) 


nnflnn 
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32=  ((2.* A*C)*B1-C*SQBT  (F))/(4.*A) 

ANUN=ABS (2. *SQFT (A) *E1*2.*A*C) 

CBN=ABS(2.*SCBT (A*F) *C) 

IF  (ANUH.LE.  1.E-04)  GO  TO  10 
IF  (DEN.LE. 1. E-04)GO  TO  10 

B3=ABS  (  (2.*SQRT  (A)  *Bl  +  2.  *A *C) /  (2.*SQBT  (A*F)  *C)  ) 

AB3=ALOG 1B3) 

AJ4=B2  +  {4«*A*F-C**2)  *AB3/ (8.*SQBT  (h**3)) 

GO  TO  11 

10  AJ4=B2 

11  B4=SQRT(E+D*F) 

B5  =  ((2.*E*D)  *B4-D*SQRT  (F)  )/  (4.*B) 

AN UH=ABS  (2.«-SQBT  (B)  *E4*2.*B*D) 

DBN=ABS ( 2. *SQRT (B*F) *D) 

IF  (ANUM.LE. 1. E-04) GO  TO  15 
IF (DEN. LE. 1.-04) GO  TO  15 

B6=ABS  (  (2. * SQRT  { B)  *  E4*2.  *B*C)/(2.*SCRT(B*F)  *E)  ) 

AB6= ALOG  (B6) 

AJ2=B5*(4.*B*F-E**2) *AB6/ (8.*SCBT (B**3|) 

GO  TO  16 

15  AJ2=B5 

16  CONTINUE 
EOT=POT1/  (2.  *AREA) 

AR  1=2.  *B*  ( AJ  1-  AJ  2)  -E*  (AJ3-AJ4) 

AB2=(2.*AH1-  (2.*  B*C—  E*D)  *  POT)  /  (4.  *  A*E-  E**2) 

POTXSI=2.*AREA*AR2 

AR3=4.*A*  (AJ3-AJ4)-2.*E*  ( AJ 1- AJ2) - (2.*A*D-E*C) 
POTETA=(2.*AREA*AR3)/(4.*A*B-B**2) 

EETORN 
END 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

CUBDIBiTC  CALCOLATE  THE  NCRRAI  VECTCB  TO  THE  SOBFACE 
,BY  LISTING  THE  EDGES  ASSOCIATED  NITB  EACH  TBIANGLE 
IN  A  SEQUENTIAL  BANNER. 

INPUT: NCONN.NBGUND, NFACE,NEEGE,NSE 
OUTPUT:  N BOUND, ININ 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  CUBDIR  (NCO  NN,NBCUND,NFACE,KEDGE,IEIN,NSE) 

c 

INTEGER  NCONN (NEDGE,3) , NBO UND  (  50, 4) ,IHIN (NEDGE) 

INTEGER  IHAX  (6 ) 

IH=  1 

IHIN(IIi)  =1 
M  1=0 

DO  999  I JK=1 ,NF ACE 
IK=I JK 
DO  2  1=1,6 
IH AX (I)  =  0 
2  CONTINUE 

IFLAG=0 
DO  4  J=  1  ,IN 

IF  (IJK  .EQ.  IMIN  ( J)  )  GO  TC  1 
IFIAG=1 

4  CONTINUE 

IF  (IFLAG.EC- D GO  TO  999 
1  J1=Q 

11  =  0 
L  1=0 

I2=NBOUNC  (IK, 2) 

I3=NBOUNC(IK,3) 
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I4=NB0UND(IK,4) 

N1=NC0NN <12, 2) 

H2=NC0NN  (12,3) 

N3  =  NCONN (13,2) 

IPf H3. EQ.N1  .OB.N3.EQ.N2)GO  TO  5 
GO  TO  10 

5  N3=NCONN  (13,3) 

10  CONTINUE 
11=1 

ITEHP=I2 

11  DO  20  IJ=1,NFACE 
DO  12  J=1,Ifl 

IP  (IJ  « EQ. ININ  ( J) ) GC  10  20 

12  CONTINUE 
J2=NBOUNC(IJ,2) 

J3=NBOUNE (IJ ,3) 

J4=NBOUND (IJ , 4) 

IP  (ITEHP.  EQ.J2  .03.  ITEHP.  EC- J3  .CB.  ITEHP.EC>  J4) 

1  GO  TO  15 
GO  TO  20 
15  IL=IJ 

GO  TO  25 

20  CONTINUE 

IP  (II  .  EQ.  1  .AND.  J1.EC.  1)60  TO  21 
J1=1 

ITPHP=I3 
GO  TO  11 

21  IP  (I1.E0. 1  .AND.  J1-EQ.1  .AND.  L1.EC-1)  GO  TO  23 
11=1 

ITEHP=I4 
GO  TO  11 

23  IP  (Hi. EQ. 1 )  GO  TO  999 

H1=H  1*  1 
IK=I JR 
GO  TO  1 

25  KN1=NC0NN  (ITEHP.2) 

KN2=  NCONN(ITEHP,3) 

IP  (N 1.  EQ. R N1  .OB.  N1.EC.KN2)  GO  TO  35 
NN3=N 1 
GO  TO  40 

35  IP (M2. EQ.KN1  .OB.  N2.EQ.KN2)  GO  TO  36 

RM3=N2 

GO  TO  40 

36  RN3=N3 

40  J2=NBOUNC(IL,2) 

J3  =  NBOUNC  (11,3) 

J4  =  NBOUN  C (IL ,4) 

IP  (J2. E0> ITEHP)  GO  TO  59 

IP  (NCONN  (J2,2).EQ.  KN1  .OB.  NCONN  (J2, 2)  .  EQ.  KN2)  GC  TO  57 
RN4=NCONK  (J2,2) 

GO  TO  68 

57  RN 4= NCONN  (J2 ,3) 

GO  TO  68 

59  IP (NCONN  (J3,2) • EQ. KN1  .OB.  NCONN (J3,2) . EQ. RN2)  GO  TO  61 

RN 4=  NCONN  (J3, 2) 

GO  TO  68 

61  KN 4=  NCONN  (J3 , 3) 

68  CONTINUE 

IPCIH.EQ.  1)  GO  TO  115 
IP (ITEHP. EQ. IflAX  (6) )  GC  TC  109 


in  A  X  (1)  =  IMAX  (5) 

IHAX  (2) =IHAX  (4) 

IH AX  (3)  =  IH  AX  (6) 

GO  TO  115 

109  IHAX(1)=IHAX  (4) 

IN  AX  (2)  =IHAX  (6) 

MAX  (3)  =MAX  (5) 

115  IP (fl  1. NE. 1 )  GO  TO  175 

IF (IT EM P. EC*  NECUND(IJK,4)  )  GO  TO  165 
IHAX(1)=NB00ND(IJK,4) 

MAX(2)  =  NBCUNE(IJK,2) 

IH  AX  (3)  UNBOUND  (IJK,  3) 

H1=0 

GO  TO  175 

165  MAX(1)  =  NBOUN£(IJK,3) 

MAX  (2)  =NEOUNE  (I JK, 4) 

MAX  (3)  =  NBCU  NE  (IJK,  2) 

El=0 

175  KD0HHY=KN3 

BO  100  1=1,2 

IF  (I.  EQ.  1  -AND.  IN.  NE.1JGC  TO  99 
ID=I  +  (1-1)  *  2 

IF  (ITENP.EQ. 12)  GO  TO  79 

IF(ITEMP.EQ.I3)  GO  TC  89 

IF  (Nl.  EQ.KN3  .ANE.  N2.  EC. Kill)  GO  TO  69 

IF (Nl.  EQ. KN1  -AND.  N2.EC.KN3)  GO  TO  69 

MAX  (ID)  =13 

IH  AX  (ID*  2) =12 

GO  TO  99 

69  IH  AX  (ID)  =12 

IHAX  (ID+2) =13 
GO  TO  99 

79  NN1  =  NCONN  (13,2) 

NN2=NC0NN(I3, 3) 

IF  (NN1.EC.KN1  .AND.  NN2.EC.KN3)  GO  TO  81 
IF(NN1.EC-KN3  .AND.  NN2-EC-KN1)  GO  TO  81 
IHAX  (ID)  =14 
IH  AX  (ID+2)  =13 
GO  TO  99 

81  IH  AX (ID) = 1 3 

MAX  (ID+2)=I4 
GO  TO  99 

89  IF (N  1.  EQ.KN3  . AND. N2. EC* KN1)  GO  TO  91 

IF  (N1.EQ.KM  .AND.  N2.EC.KN3)  GO  TO  91 
IHAX  (ID)  =14 
IH  AX (ID+2) =12 
GO  TO  99 

91  IN  AX  (ID) =12 

MAX  (ID+2)  =14 

99  KN  3=KN4 
12= J2 
I3=J  3 
14  =  J4 

100  CONTINUE 
KN3=KDUN  BY 

NA  1=NCONN  (MAX  (1)  ,2) 

NA 2=  NCON N (IHAX  (1) ,3) 

SB1  =  NC0NN  (MAX  (4)  ,2) 

NB2=NCONN  (IHAX  (4) ,3) 

IF (NB 1. EQ. NA  1  . OR. NB 1 . EQ . NA 2)  GO  TO  125 
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IF  (HB2.  EQ. NA 1  .OB.  NE2.E£.NA2)  GO  TO  125 
IDUHHI=IHAX  (6) 

IHAX  (6)  =  IHAX  (4) 

IN AX  (4) =IDUHH Y 
125  INAX  (2)  =ITEMP 
IHAX (5) =ITENE 
IF(Ifl.NE.I)  GO  TO  149 
NBODND(IK,2)=IHAX  (1) 

NBOONC  (IK,  3)  =  IHAX  (2) 

NBOONO  (IK  ,4)  =IMAX  (3) 

149  NBOUND(Il,2)=IHAX(6) 

NBOUND(IL,3)  =  MAX  (5) 

N'BOUN  D  ( IL,  4)  =  IHAX  (4) 

IH=IM+1 
IHIN(IH) =IL 
IK=It 

IF  (M.  EQ.  NFACE)  GO  TO  1000 
GO  TO  1 

999  CONTINUE 

1000  CONTINUE 
IF(NSE.EC.O)  GO  TO  1001 
SBITE<21,  98) 

98  FOFMAT  ( *  1*  ) 

NRITE  (21  ,  102) 

102  FORMAT  (10X, ‘LIST  OF  EDGES  £  VERTICES  HOUNDING  EACH 

1  FACE*) 

DO  1999  IJK-  1 , NFACE 
I2  =  NBOUNE(IJK,2) 

I3=NB0UNE(IJK, 3) 

14 UNBOUND (IJK #4) 

IF  (NCONN (12, 2) . EQ. NCONN (13,2))  GC  TO  1005 
IF  (NCONN  (12,2) .EQ. NCCNN (13,3) )  GO  TO  1005 
Nl^NCONN  (12,3) 

GO  TO  1006 

1005  N1=NC0NN(I2,2) 

1006  IF (NCONN <13, 2) .EQ. NCONN  (14,2) )  GC  TO  1010 
I? ( NCONN (13,2) . EQ. NCONN (14,3) )  GO  TC  1010 
N2S NCONN  <13, 3) 

GO  TO  1011 

1010  N2S NCONN  (13,2) 

1011  IF  ( NCONN  (14,2)  .EQ.  NCONN  (12,2) )  GO  TC  1015 
IF  (NCONN  (I4,2).EQ. NCONN  (12,3))  GO  TO  1015 
M3=NCONN(I4,  3) 

GO  TO  1016 

1015  N3sNCONN  (14, 2) 

1016  CONTINUE 

WRITE  (21,  10  50) IJK,I2,I3,I4,N1,N2,N3 
1050  FOHHAT(/3X,*FACE*  ,13, IX,  »IS  BOUNDED  BY  EDGES*  *  IX, 

1  13, IX, 13, 1X,I3,2X,'AND  VERTICES* , 1 X ,13 , 1 X, 13 , IX, 13) 

1999  CONTINUE 

1001  RETURN 
END 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

INTGRL:CALCULATE  1/R  INTEGRATION  OVER  A  SOURCE  TRIANGLE 
INPUTS:  (XS  (3)  ,ZS  (3)  )  :CCORS.  OF  NODES  OF  THE  TRIANGLE 
X  F,Z  F: CCCR .  OF  FIELD  POINT 
OUTPUTS: EOT=  FOT1 :TNTEGAL  RESULT 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  I NTGBL  (XS , ZS , X E , ZF , POT) 

DIMENSION  XS  (3)  ,ZS  (3) 

COMMCN/POTEN/POT 1 

FIND  THE  ENDS  OF  EACH  E  A  EG  E  ( P  ROM  (XM,ZM)  TO  (XP,ZP) 

THE  EDGE  NO  IS  DEFINED  AS  1,2,3  BETWEEN  NODES  ( 1 , 2)  ,  (2 , 3)  #  (3,  1)  A 

AND  THE  DIRECTION  IS  DEFINED 

AS  FROM  NODES  1- 2,2-3, 3- 1, RESPECTIVELY. 

POT=0. 

DO  10  1=1,3 
XM  =  XS  (I) 

ZM=ZS  (I) 

11=1*1 

IF  (1 1  .EC.  4)  11  =  1 
XP=XS (II) 

ZP=ZS  (II) 

XPXF=XP-XF 
ZPZF=ZP-2F 
XMXF=XM-XF 
ZHZF=ZM-ZF 
XPIM=XP-XM 
ZPZH— ZP—  Z  M 

AR=ABS (XPX?*ZMZF-XMXF*ZP2F) 

IF  (AR  .LE.  1.E-12)  GO  TO  10 
XPXM2=XPXM**2 
ZPZM2=ZPZM**2 
EL=SQRT  (XPXM2*ZPZfl2) 

F0=AR/DL 

EL2=DL**2 

XPXHIL=XPX112/CL2 

ZPZHL=ZPZM2/EL2 

DOT=  (-XfiXF*XPXH-ZHZF*ZFZE)/EL2 

XO=Xfl*DOT*XPXM 

ZO=ZM*DOT*ZFZK 

P02=R0**2 

DLP=  ((XP-XO)  *XPXfl*(ZP-ZO)*ZFZn)/EL 
DLI!=  ((XM-XO)  *XPXM+(ZM-Z0)*ZEZ!1)/DL 
BP=SQRT  (E02*ELP**2) 

Rfl  =  SQRT(R02  +  DLM**2) 

SIGN=-  ((XO-Xi)  *ZPZM-(Z0-ZF)  *XPXf!)  /  (  CL*F0) 

RATIO=  (RP*DLP)/(HM*DLM) 

FOT=POT*SIGN*FO*  A LOG  (RATIO) 

0  CONTINUE 

FOT=ABS  (EOT) 

EOT  1  =  POT 
EETURN 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  MATBIX:TO  OBTAIN  SOURCE  VECTORS  CIA, CIB, COUPLING  HATRICES  T,TT. 

C  THE  WIRE  IS  INFINITELY  LONG  OR  LOADED.  THE  APERTURE  IS  OF 

C  ARBITRARY  SIZE  AND  SHAPE, AND  IS  IN  AN  INFINITE 

C  CONDUCTING  PLANE  OF  ZERO  THICKNESS. 

C  INPUT: NA,NE:TOTAL  *  OF  EXPANSIONS  IN  APERTURE  AND  OH  HIRE 

C  NNODE, NEDGE, NPACE:  *  OF  NODES, EDGES , FACES  OF  APERTURE 

C  DATNOB  (NNODE,3)  :NOCE  I,  X,Z,  CCOB.  OF  ALL  NODES 

C  NCONN  {NEDGE,  3)  :  EDGE  #,  NODE  #(FBOH,TO)OF  ALL  EDGES 

C  NBOUND  (SO, 4) : FACE  t , E DGE  «  OF  ALL  TRIANGLES 

C  ETHETA,EPHI, THETA, PHI: INCIADENT  E-FIELD  ARPLITUDES  6  ANGLES 

C  OUTPUT:  CIA  (N  A)  ,CIB  (NA)  ,T  (N  A,  NB)  ,  TT  { NB,N  A)  ,  Y  (NA,NA) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

SUBROUTINE  H AT RIX  (N A, NE, KNOB E, NEDGE, N FACE, CAT NOD, NCONN, 

1  NBOUND, CIA, CIB, T,TT,Y) 

COMPLEX  Y  ( NA  , N A)  ,CIA  (NA)  ,T  (NA,NE)  ,  HPHI,  HTHETA 
COMPLEX  H1,HX,HZ,FX, FZ,TF, SPOT,CPHI , JK,CAXSI ,C A ET A 
COMPLEX  CIB(NA) ,TT  (NB, NA) 

COMPLEX  C1,C2,C3,C4,C5,G1,G2, V 1 , V2 , AIP , AIM ,TP, TH,T 1 (19) ,TNB(19) 
REAL  DATNOB  (NNODE, 3)  ,XF(3)  ,ZF(3)  ,LF(3)  ,XS(3)  ,ZS(3)  ,LS(3)  ,CT(3,2) 
INTEGER  NCONN (NEDGE, 3) , NBC UND (50,4) 

INTEGER  £F(3),NF  (3)  ,  ES  (3)  ,NS  (3) 

COHHON/L HIRE/C, HL 

COHMON/FIELD/HPHI, HTHETA, PHI, THETA 
COMMON/KKK/AK , PI 
COMMON/JKK/JK 
COHMON/H  EU/H  E, HU 
COMHON/LCAD/G 1 ,G2 
COMMON/VOLT/ V 1, V  2 
COMMON/ZOO/ZO 
C 

C  NSE=TOTAL  NO.  OF  SURFACE  EDGE  (BOUN DAB Y)  ;  NA-TOTAL  NO.  OF 

C  INTERNAL  EDGES 

C 

NS E=NEDGE-NA 
DO  60  M=  1 ,  NA 
BO  70  N=  1, NA 
70  Y(H,N)  =  (0.,0.) 

Cl  A  (M)  -  (  0.  ,  0.  ) 

CIB  (M)  =  (0. ,0.) 

T 1 (M) =(0.,0.) 

TNB (M)  =  (0. ,0. )  N 

DO  80  N=  1, NE 
IT(N,M)  =  (0.  ,0.) 

80  T(M,N)  =  (C.,0.) 

60  CONTINUE 

COST=COS  (TFETA) 

SINT=SIN  (THETA) 

COSP=COS (PHI) 

SI NP=SIN  (PHI) 

HX=HTHETA*COST*COSP-HPHI*SINP 

HZ=-HTHETA*SINT 

C4*G  1*G2 

C5a  1.-C4 

C 1 =~G  1/C  5 

C2=C4/C5 

C3=-G2/C5 

C4*2. *C5 

AIPa((1.-Gl)  •VUG1*(1.-G2)*V2)/C4 


AIMS(-G2*(1.-G1)*V1-(1.-G2)*V2)/C4 


9*? 


V? 


non  on  onnoo  noon  nnonnno  nonooo 


DO  10  I  J=  1 , NF  ACE 
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FOB  EACH  FIELD  TBIANGLE  IJ, OBTAIN  TBIANGLE  PABABETEES: 

EDGE  #:  EF (3)  ;  NODE  *  :NF(3);  X.Z.  COOB.  OF  MODES  XF(3),ZF(3); 
LENGTH  OF  EDGES:  LF(  3) 

OBTAIN:CENTBOD:XC,ZC;  TESTING  COBFONENTS  OF  EDGES: CT (3,2) 

N1 1=IJ 

CALL  TPAEB  (N  1 1 , D ATNCD, NCON N, NBCUND , NNCDE , NEDGE , 

1  EF,NF,XF,ZF,IF) 

XC=  (XF ( 1 ) +XF  (2) +  XF<3) ) /3. 

ZC=  (ZF  (1)  *ZF  (2)  *ZF(3)  )/3. 

CT  (1,1)=  (XF(2)  ♦XF(3)  )/2.-XC 
CT  (1,2)=  (ZF  (2)  *ZF(3)  )/2.-ZC 
CT  (2,  1)  =  (XF  ( 3)  ♦XF  ( 1) )  /2«-XC 
CT  (2,2)=  (ZF  (3)  *ZF  (1)  J/2.-ZC 
CT  (3,  1)=(XF(1)  *XF  (2) )  /2.-XC 
CT  (3,2)=  (ZF  (1)  *Z  F  (2)  )  /2. -ZC 
DO  40  IJK= 1, NFACE 

FOB  EACH  SO 0 BCE  TBIANGLE  IJK: 

OBTAIN:  EDGE  PABABETEES :  E  CG  E  *  ES(3),NCDE  *  NS(3),COOB.  OF  NODES 
XS(3),ZS(3);  LENGTH  CF  EDGES  IS  (3) 

OBTAIN: ABEA  (ABEA  GE  TBIANGLE  IJK) 

BT  TAKING  MAGNITUDE  OF  VECTCB  CBCSS  PEODDCT  CF  TNG  SIDES 
K1  1=IJK 

CALL  TPAEB  (N  11, DAT  NOE,  NCCNN,  NBOUND 
1  ,NNODE,NFDGE,ES,NS,XS,ZS,IS) 

ABEA=  (XS  (2)  -*XS  (1))*  (ZS  (3) -ZS  (  1))  -  (XS  (3) -XS  ( 1)  )  *  (ZS  (2)  -ZS|1)) 

AB  EA=  ABS  (AFEA)/2. 

OBTAIN  SCALAfi  &  VECTCB  POTENTIAL  INTEBALS ( 1/2* ABEA  EXCLODED)  : 

CP HI  ,  CAXSI ,CAETA 

CALL  SCAINT (XS,ZS,XC,ZC,CPHI, AREA) 

CALL  VECINT  (XS, ZS ,XC, ZC, CAXSI, CAST A , ABEA) 

DO  20  IB=1 ,3 

FOB  EACH  MODE  IR  OF  THE  FIELD  TBIANGLE  IJ: 

OBTAIN: FIELD-EDGE-TESTING  INDEXsHjDIBJECTICN  OF  COBBEMT: DIBJ 


N1 1=IB 

CALL  TADJ  (N1 1,NA, EF, NF, NCC NN, NEDGE, B ,DI&J) 

IF (DIBJ  .EC.O.)  GO  TO  20 
IF  (I«JK.  NE.  1)  GO  TO  1 

COMPUTE  SOUBCE  VECTOR  CIA(R):  INCIDENT  PLANE  HAVE  IN  REGION  A 
BC=XC*SINT*COSP*ZC*CCST 

HI  =  (CT  (IR,  1)  *HX*CT (IR,2)*HZ)*CEXP { JK*8C) 

CIA  (B) =CIA  (M) *2. (IB) *H1 

COMPUTE  COUPLING  MATRIX  T(fl,N),  SOURCE  CIB(B)  FOR  TEH  IN  REGION  B 

E22=D**‘2*XC**2 

TL=HL 

lL2*TL/2. 

ZCT*ZC»TL2 

B1=SCHT(ZCT**2+E22) 


3^ 

Hi  =  (  1.  -ZCT/R  1)  *C  EXP  (-JK*  (B  1*TL2)  ) 

D4=LF(IR)*DIRJ*B*CT (IP, 1)/ (2.*PI*C22) 

TP=D4*2. *CEXP ( JK*ZC) 

TM=D4*2.  *CEXP  (— JK*ZC) 

11  (M)=C1*TM*C2*TP*T1  (M) 

TNB  ( M) =C2*TM ♦C3*TP*TNB (M) 

T ( H,  1) =T  |H,1)*D4*H1 

ZCT=ZC-TI2 

R1=SQBT (ZCT**2»D22) 

H 1  =  (  1.+ZCT/R  1)  *CEXP(-JK*  (RUTL2)  ) 

T(M,NB)  =1  (M,NE)  ♦  C4* H 1 
CO  30  N 1  =  2, NE-  1 

ZN=  WL*FLCAT ( 2*  N 1- N8- 1 ) /E LC AT (2*NE-4) 

BN  =SQRT  (  C22*  (ZC-ZN)  **2) 

Hl=(JK/RN**2*1./RN**3)  *CEXP  (- JK*  BN)  *HL/FLO AT  (NB-2) 

30  T(M,N1)=T(M,N1)  *D4*D22*H1 

CIB  (M)=AIP+TM*AIM*TP+CIE (B) 

C 

C  COMPOTE  ADMITTANCE  MATRIX  Y  (M, N) 

C 

1  CONTINUE 

CO  50  IK=  1, 3 
C 

C  FOB  EACH  NODE  IK  OF  THE  SOURCE  TBIANGLE  IJK: 

C  OBTAINiSOUBC E-EDGE- INDEX  Nj  CURRENT  DIRECTION  DIRS 

C 

Nl  1=IK 

CALL  TADJ (Nl 1,NA, ES, NS, NCONN, NEDGE, K ,EIBS) 

IF  (DIRS  -EC-0-)  GO  TO  50 
C 

C  COMPUTE  SCALR  POTENTIAL  SPOT, VECTOR  POTENTIAL  IN  X.  Z.  :FX,PZ 

C  ;T HE  DOT  PRODUCT  OF  P01EN1IALS  6  TESTING, IF, IS 

C 

A2=DIRS*LS (IK)/(4.*PI*AREA*2.) 

FX=A2*((XS  (11-XS  (IK))*CPHI*(XS(2)-XS(1))*CAXSI 
1  ♦  (XS(3)-XS(1|)*CAETA) 

FZ-A2*  ( (ZS  (1)-ZS  (IK)  )  ♦CPHI*  (ZS  (2)-ZS  (1)  )  *CAXSI 
1  *(ZS(3)-ZS(1))*CAETA) 

TF=FX*CT  (IR,  1)  *FZ*CT  (IB, 2) 

SPOT=— CPHI#A2*2.  /(0.,1.)/WU 

Y ( M, N)  = Y (M, N) ♦4.*DIHJ*LF (IB) *  (  (0. , 1. ) *NE*TF-SPOT) 

50  CONTINUE 

20  CONTINUE 

40  CONTINUE 

10  CONTINUE 

C 

C  GALERKIN  SOLUTION  :Y(H,N)=Y(N,M) 

C 

CO  90  H= 1, NA 
DO  90  N=  1 ,  M 

Y(  H,  N)  =  (Y  (M,  N)  ♦  Y(N,M))/2. 

90  Y(N,M)=Y  (M,N) 

C 

C  CBTAIN  TT  (M, N) 

DO  100  N  =  1 , N  A 
CO  100  M=1,NE 
100  TT  (M,N)=-T  (N,H) 

IF  (G 1  .EQ.  0.  .AND.  G2  .EC.  0.)  GO  TO  2 
C 

C  ADD  T1 ,TNB  DUE  TO  LOADS  TC  MATRIX  T: 


c 
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CO  110  M= 1 #  N A 
T(H,1)=T(H,1)*T1  (N) 

110  T(fl#N3)=T  (fl,  NE)  *TNB  (H) 

2  CONTINUE 

BETUPN 
END 


J*fm  "*•  •' 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  It!  P:  FIN  D  IMPEDANCE  MATRIX  Z  FOR  AN  I N FI  NIT ELY-LONG  WIRE 

C  INPUT:NB:#  EXPANSIONS. 

C  RB: RADIUS  OF  WIRE 

C  OUTPOTlZ 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  IMP(NE,RE,Z) 

COMPLEX  Z(N6,NB)  ,ZHN(5)  , P , JK, E, El , E2 
DIMENSION  CZ  (5) 

COMHON/K  KK/AK , PI 
COMMON/JKK/JK 
COMMON/LWIRE/D, WL 
COHMON/W  EU/WE ,  WU 
*LB=WL/FLOAT(2*(NB-2)  ) 

C 

C  FIND  Z(M,N),  M,N={2,  NE-1) 

C  H=2  ,  THE  1ST  SUESECTICN  CN  WL,  SO  M-1  IN  MO,N-1  IN  NO 

C  (MM,  MO,  MP)  ,  (NM,  NO,  NP)  ARE  THE  POINTS  ON  SUBSECTION  MSN 

C 

C2  =  D*2. 

NM=0 

NO=1 

NP=2 

H1  =  2 

M2  =  NB-1 

CO  10  M  =  M  1 ,  N  2 

HO=2*  (M-l)-l 

HP=MO>1 

HM=MO- 1 

CZ  (  1)  =ABS (FLOAT (MO-NO)  *ALB) 

DZ  (2)=ABS(FLCAT(MP-NP)*AIE) 

CZ  (3)=ABS(FLOAT(MP— NM)*AIB) 

DZ  (4)=ABS(FLCAT(HM-NP)*AIE) 

DZ  (5)  =  ABS  (FLOAT  (MM-NM)  *AIB) 

DO  20  J=  1 ,5 

20  ZMN(J)=P  ( ALB , DZ ( J)  ,R B) - P ( ALB , DZ { J)  ,D2) 

Z(H,M1)=(0«,1»)*NU*4*A1B**2*ZMN(1) 

Z(H,Ml)— Z(M,  Ml)-(0.,1.)/NE*(ZMN(2)— ZHN{3)— ZMN(4)»ZMBf5)) 

Z ( M 1  ,M) =Z  (M, H 1 ) 

DO  30  J=  1,  M2-M 
Z<H*J,H1+J)=Z  (M , H 1 ) 

30  Z(H1*»T,M*J)=Z(M,M1) 

10  CONTINUE 

C 

C  FIND  Z(M,1)  =  Z(1,M)  ,Z  (M,NE)=Z(NB,M)  ,  BY  USING  PULSE  TESTING 

C  OF  J  (M)  ,  M=  (  2,  NB—  1) 

C 

ALS-ALB 
E22=D2**2 
L=  (NB-D/2+1 
1L  =  WL 
TL2=TL/2. 

E=CEXP(-JK*TI2) 

CO  40  M=2,L 

Z0=W1*FLCAT (2*M-NB-1) /FLOAT  (2*NB-4) 


CZZ={ZO) 

DZ 1sABS  (CZZ+TL2) 

DZ  2s  ABS  (  CZZ— TL2) 

ZH»(1)  =  P  (ALB  ,DZ  1,RB)-P(ALB,CZ1,D2) 
ZMN  ( 3)  =P  ( ALB,DZ2,RB)  -  P(ALB,DZ2,D2) 


'j  -■ 


n  n  n 


ZS=  (ZO-ALS) 

ZPl=ABS(ZP+T12) 

ZP2=ABS  (ZP-TI2) 

ZN1=ABS(ZM+TI2) 

ZM2=ABS  (ZM-TI2) 

ZNN  (2)  =  P  (ALS  ,ZP1,RB)-P  (AIS,ZP1,D2) 

1  -P(ALS,ZM1,EB)*P(ALS,ZN1,D2) 

ZHN(4)=P{ALS,ZP2,RB)— P (AIS,ZP2,D2) 

1  -P  (ALS,Z(12,BE)  ♦P(ALS.Z«2,E2) 

E1=AK/HE*2.*AIB*E 
E2=-(0.#  1.)/KE*E 
Z(H,  1)  =E  1*ZHN  ( 1)  ♦E2*ZMN  (2) 

Z<M,NB)=£1*ZMN  (3) -E2*ZMN  (4) 

Z(1.M)=Z<f1,1) 

Z  ( NB  ,  fl)  =Z  (M,  NE) 

K1=NB-M+  1 
Z(M1,1)  =  Z(M,NB) 

Z(fl1,NB)  =Z  (il*  1) 

Z(1,M1)  =  Z(M1,1) 

Z(NB,M1)  =Z  (M  1.  NB) 

40  CONTINUE 

FIND  Z  ( 1  ,  1 )  =  Z  (  NE  r  NB)  ,  Z  ( 1#NE)  =Z  (NB,  1) 
TL22=TL**2 

U01=TL*SCBT(TL22*RB**2) 
002=TL+SCBT(TI22>D22) 

CALL  SICI{S1,C1,U01*AK) 

CALL  SICI(S2,C2,U02*AK) 

E=CEXP{— JK*TL) 

E1=E/HE 
I2=(0.,  1.)  *E1 

Z(1.1)=-E2*(E(ALS,0.,BB)-P(ALS,0.,C2)) 

ZC  1#NB)  =  E2*<P(ALS,TL#BB)-P(AIS,TL,D2|) 

1  ♦1./ME*AK/(2.*PI)*<-C1*C2*  (0.  ,  1.  )  *  (£  1-S2) ) 

Z{NB,NB)  =Z(1,  1) 

Z ( NB ,  1 )  -Z  { 1  #  KB) 

BETURN 
END 


S' 

>v> 

& 


r-> 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  SICIiCOMPOTES  SINE  and  COSINE  INTEGRALS.  INPOT  IS  X; 

C  OOTPOTS  ARE  SI,  Cl. 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  SICI { SI, Cl , X) 

Z=ABS  (X) 

IP  (Z-4.  )  1,  1,  4 
1  1=  (4.-Z)  ♦  (4.+Z) 

SI  =  X* ( (( Cf1.  753141E—  9*Y  +  1.568988E-7) *Y+ 1. 374168B-5) 

1  *Y+6.939889E-4) *Y+ 1. 964882 E-2) *Y+4. 395509E-1) 

CI=(  (5.77215EE-1+ALOG (Z») /Z-Z*  (  ( ( ( { 1. 386985B- 10*1 

1  *1. 584996 E-8) *Y* 1-7257 52 E-6) ♦Y+1. 185999E-4) *Y+ 

2  4.  990920E-3)  *Y+1. 315308E-  1)  )  *Z 
BETURN 

41  SI=SIN  (Z ) 

Y=  COS (Z) 

Z=4./Z 

0=  (((((( {(4.  048069E-3*Z-2.279143E-2) *Z+5. 515070B-2) 

1  *Z-7.  261642E-2)  ♦Z«-4.9877  16E-2)  *Z-3. 3325 19E-3)  *Z- 

2  2.  3 146 17  E- 2)  *Z-  1.  1  34958E-5)  *Z*6.2500 1 1E-2)  *Z+ 

3  2.583989E-10 

7=  <((((( <(<- 5.  108699E-3*Z+2.819179E-2)*Z— 6. 537283E— 2) 

1  *Z  +  7.902034E-2) *Z-4.400416E-2) *Z-7.945556E-3) *Z+ 

2  2- 60 1 293 E- 2) *Z- 3. 7640  COE-4) *Z-3. 1224 18E-2) *Z- 

3  6. 646441 E- 7) ♦Z+2.500000E-1 
CI=Z*  (SI*V-Y*U) 

SI=-Z*(SI*U  +  Y*V)  +1.570796 
BETORN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  Pi  FIN 0  THE  INTEGRAL  OE  1/  (4*PI*2*AL)  *CEXP  (-JMH)  /B  IN  A 

C  SUBSECTION  OF  S*AL 

C  Z— DI ST.  IN  Z  BETWEEN  FIELD  6  SOUBCE  PTS. 

C  ZL-DIST.  IN  ZL (=SQRT  (X**2+Y**2) ) 

C  AL=H ALF  SUBSECTION  LENGTH 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
COMPLEX  FUNCTION  P{AL,Z,ZL) 

COMPLEX  PO , JK 
BEAL  11,12,13, 14, LR,KL 
COMMON/KRX/AK, PI 
COHMON/J  KK/JK 


C 

c 

c 


fi-SQBT  (ZL**2+Z**2) 

F0=CEXP(-JK*B)/(4.*PI) 

IP (R  .LT.  (  10* AL) )  GO  TO  30 

FOB  R  >=  10*AL,  FIND  P=P  ( AO, A 1 , A2 , A3 , A4) 

KL=AK*AL 

1R=AL/B 

ZR2= (Z/R) **2 

ZR4=  (Z/R ) **4 

100=-  1.  +3.  *ZR  2 

A01=3.-30.*ZR2+35.*ZB4 

A0®1.+LR**2/6.*AOO+LR**4/40.*A01 

A1®LR/6.*A00+LR**  3/40.  MCI 

A2=-ZR2/€.-LR**2/40.*(1.-12.*ZR2+15.*ZB4) 

A3*LB/60.*  (3.*ZR2-5.*ZR4) 

A4=ZR4/120. 

P*P0/R* (A0+(0.,1.)*KLM1+KL**2M2  +  (0.,1.) *KL**3M3*KL**4M«) 
GO  TO  40 


39 


FOB  R<  10*AL.  FIND  P=P  (11,12,13,14)  : 


30  Z1=Z*AL 

Z2=Z-AL 

G1=SQRT (ZL**2*Z 1**2) 

G2=S0RT (ZL**2*Z2**2) 

C 

C  FOR  Z<=  AL: 

C 

I1=AL0G(  (Zl+Gl)*  (-Z2*G2)/ZI«'*2) 

IP (Z  .LE.  AL)  GO  TO  5C 
C 

C  FOB  Z  >  AL: 

C 

I1=ALOG(  (Z1*G1)/(Z2*G2) ) 

50  12=2. *AL 

I3=Z  1/2.  *G1-Z2/2.*G2*ZL**2/2.*I1 
I4=2.*AL*ZL**2*  (2.*AL**3*6.*AL*Z**2)/3. 

P=P0/(2.*AL)  *(I1-JK*(I2-R*I1) 

1  -AK**2/2.*(I3-2.*R*I2*R**2*I1) 

2  *(0., 1.) *AK**3/6.*  (I4-3.*R*I3*3.*B**2*I2-R**3*I1)) 

40  RETURN 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  CSHINVrlNVERTS  MATRIX  A  Of  DIMENSION  NXN  (N=BDIH),ANC  STORES 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  BESULTS  IN  A 

SUBROUTINE  CSMINV (A,NDIM,N) 

COMPLEX  A ( NDIrt, NDIM) , PIVOT (60) , AM AX, T,SUAP, DETERM 

COMPLEX  U,CMPLX ,CON JG 

INTEGER*4  IPIVOT  (60)  .INDEX  (60,  2) 

BEAL  TEMP, ALPHA (60) , CABS 
COMPLEX  CTEME, CALPHA (60) 

IEER=0 

IP (NDIM  .LE.  60)  GO  TO  5 
IERR= 1 

hRITE  (3 ,4)  NDIM 

4  FORMAT  (*OCSMINV  ERROR,  ATTEMPT  TO  INVERT  A  MATBIX* 

1  14.*  ON  A  SICE, * / '  WHEN  60  X  60  IS  THE  MAXIMUM  ALLOHED*) 

RETURN 

5  CONTINUE 
DBTERM=CMPLX (1.0,0. 0) 

SUHAXA=0. 

CO  20  .7=  1,  N 
ALPHA  (J)=0.0 
CALPHA  (J)=  (0.0, 0.  0) 

SUHRON=0. 

CO  10  1=1,  N 

CALPHA (J) =CALPHA (J) ♦A (J,I) *CONJG (A(J,I)  ) 

ALPHA  ( J ) =REAL (CALPHA (J)  ) 

10  SUHROH=SUMROk  +  CABS(A  (J,  I)  ) 

ALPHA  (J)  =SQRT  (ALPHA  (J)  ) 

IF  (SUM RON  .GT.  SUMAXA)  S UMAX A=SUMRON 
20  IPIVOT (J)=0 

SO  600  1*1, N 
AMAX=CMPLX(0.0,0.0) 

CO  105  J=  1  ,N 

IP  (IPIVOT  (J)  -1)  60,  105,60 
60  CO  100  K= 1 , N 


IP  (IPIVOT (K) -1 )  80,100,740 
80  CTEHP= AHAX*CCN JG  {  AHAX)  -  A  (J,  K)  *CONJG  (A  (J,K)) 

TEHP=REAL  (CTEHP) 

IP  (TEMP)  85,85,  100 
85  IROH=J 

ICCL0H=K 
AH AX=A  ( J  ,  K) 

100  CONTINUE 

105  CONTINUE 

IPIVOT  (ICOLUH)  =1 PIVOT (ICCLUH)  »1 
IF  C  IROU  -  ICC1UH)  140,260,140 
140  DETERH=  -DETERH 

EO  200  L= 1 , N 
SB  AP=A (IEOH, L) 

A(IROH,L)=A(ICCIUH,L) 

200  A ( ICOLUH ,L) =SH  AP 

SVAP-ALPHA  (I ECU) 

ALPHA  (IROH)  =  ALPHA  (ICCLUH) 

CALPHA  (ICOLUH) -SWAP 

AL  PH A (ICOLUH) -REAL(CALPHA (ICOLUH)  ) 

260  INDEX  (I,  1)  =1  BOH 

INDEX (I, 2) SICCLUH 
PIVOT  (I)  =A  (ICCLUH,ICC1UH) 

0=  PIVOT  (I) 

ALPHAI= ALPHA  (ICCLUH) 

CALL  DTRHNT  (CETERH, U,ALPRA I) 

CTEHP=PI  VOT  (I)  *CCNJG  (EIVCT  (I)  ) 

TEHP=REA1 (CTEHP) 

IP  (TEMP)  330,720,  330 
330  A ( ICOL UH, ICOLUH)  =0(1  PLX  (1.0,0.) 

CO  350  L=1,N 
U=PIVOT  ( I) 

350  A (ICOLUH ,L)  =  A  (ICOIUH,L)/D 

380  EO  550  L  1=  1, N 

IF  (LI  -  ICCLUH)  400,  550,400 
400  T=A(L1, ICOLUH) 

A (LI , ICOLUH) -CHPLX (0. 0,0.0) 

EO  450  L= 1, N 
U=A  (ICOLUH, L) 

450  A(L1,L)=A(L1,L)-0*T 

550  CONTINUE 

600  CONTINUE 

62  0  DO  710  1=1,  N  1 

L=N+  1-1 

IF  (INDEX  (L,1)  -  INDEX  (1, 2) )  630,710,630 
630  JRON=INDEX  (L,  1) 

JC OLUfl=I NDEX  (L , 2) 

CO  705  K=1,N 
SB  AP=A  (K  ,  JEGH) 

A (K, JROW) =A (K, JCOLUH) 

A(K,JCOLUH)=SHAP 
705  CONTINUE 

710  CONTINUE 

SUHAXI-0. 

EO  910  1=1, N 
SOHBOH=0. 

DO  900  3=1, N 

900  SUHROH=SUHROH*CABS(A(I,J) ) 

IP  (SUHRO  K  .GT.  SUHAXI)  SUHAXI=SUHFOH 
910  CONTINUE 


RETURN  41 

720  WHITE  (3,730) 

730  FOBMATfO*  ,10  (•******•)/ »OMATRIX  IS  SINGULAR  V*  0* 

1  ,lO(********n 

740  fiETURN 

END 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  ETRMNT: 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  DTRMNT (DETERN, U,A) 

COMPLEX  CETERA, U,CNPLX 
BEAL  CABS 

COHMON  /SC A FAC/ IS CAL  E 
BATA  ISCALE/C/ 

IF  (CABS  (CETERA)  -GT.  1.E-10)  GO  TO  100 
DETERa=DETERM*1. E10 
ISCALE=ISCALE* 1 

100  BETERM=DETERM*U/CAPLX ( A, C. 0) 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  GAUSS:SOLVE  FOR  X :  A  (N,  N)  *1  ( R)  =E  (  N)  ;  STORED  IN  B(R) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  GAUSS(N,A,B,EFS,ISH) 

COMPLEX  A  (N,N)  ,E  (N)  ,C,T 
NM 1=  N-  1 
DO  10  K=1,NM1 
C=  (0.0, 0.0) 

BO  2  I=K,N 

IF  (CABS  (A  (I#  X)  )  . LE.CABS  (C)  )  GO  TO  2 
C=  A  (I,K) 

10=1 

2  CONTINUE 

IF  (CABS  (C)  -GE-EPS)  GO  TO  3 

IS»=0 

RETURN 

3  IF(IO-EQ-K)  GO  TO  6 
BO  4  J=K , N 
T=A(K,J) 

A{K,J)=A(I0,J) 

4  A(I0,J)=T 
T=B(K) 

B(K)=B(I0) 

E(I0) =T 

6  KP 1=K* 1 

C=  1./C 
E(K)=B(K)*C 
IX)  10  J=KP  1#  N 
A(K,J)=A  (K,J)*C 
CO  20  I=KP1,N 

20  A(I,J)=A|I,J)-A4I#K)*A(K,J) 

10  B(J)=B(J)-A(J,K)  ♦B(K) 

B(N)=B(N)/A(N,N) 

CO  40  K=  1,  NH  1 
I=N-K 

C=  (0-0, 0.0) 

IP  1=1*1 
CO  50  J=IP  1,  N 
50  C=C*A  (I,  J)  *8  (J) 

40  E(I)*B(I)-C 


U2 

IS  B=  1 

RETURN 

END 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  HUL:C(L,  N)  -A  (L,M)  *B  (M,  N) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  MUL(L,H,N,A,B,C) 

COMPLEX  A  (L,M)  ,E  (M,N)  ,C  (L,N)  ,W 
EO  20  1=1,  L 
CO  20  K=  1 ,  N 
»=  (0,0) 

EO  10  J=  1 , M 

10  H=A(I,J)  *B(J.K)  ♦  « 

20  C  (I ,  K)  =H 

ffETURN 
END 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  KUL1:  A(L,t1)*E  <H)=C(L) 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  MU1 1  ( L ,M, A, E , C) 

COMPLEX  A  (L,  M)  ,  B  (M)  ,C  (L) 

CO  10  1=1,  L 
C(I)  =  (0.,0.) 

CO  10  J=  1 ,  M 

10  C(I|=A(I,J)  *B{J)  *C(I) 

RETURN 

END 
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